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SUMMARY 


The  report  analyzes  by  means  of  examples  the  applicability 
of  the  finite  element  concept  to  hyperbolic  equations.  One  of 
the  aims  is  to  provide  an  appreciation  of  the  basic  phenomena. 

The  elements  chosen  are  quadrangles  with  bilinear,  biquadratic 
or  bicubic  shape  functions.  The  equation  discussed  most  is 
the  two-dimensional  wave  equation.  Variational  approaches  are 
rejected  because  they  refer  to  problems  that  are  not  well  posed; 
instead  the  method  of  weighted  residuals  with  rather  general 
weight  functions  is  adopted.  The  equation  of  convergence  which 
poses  difficulty  if  an  extremum  formulae  exists,  has  particular 
importance  for  hyperbolic  problems  because  the  exact  solutions 
are  not  damped  and  truncation  errors  may  cause  a  neutrally  stable 
solution  to  become  unstable.  Besides  staibility  the  question  of 
accuracy  is  explored.  It  determines,  for  a  stable  method,  the 
admissible  element  size.  Hyperbolic  problems  have  equal  wavi¬ 
ness  in  the  space  and  in  the  time  direction.  Therefore,  one 
has,  from  the  point  of  view  of  accuracy,  a  natural  limitation 
of  the  Courant  number  to  1,  even  if  this  limitation  is  not 
needed  for  stability  reasons.  First  a  number  of  semi-discretized 
approaches  (discretized  in  space  but  not  in  time)  are  investigated. 
The  solutions  always  have  dispersive  character  (no  damping,  but 
wave  velocities  different  from  the  actual  one) .  The  ratio  of 
the  wave  velocities  obtained  from  the  original  partial  differen¬ 
tial  by  the  numerical  approach  provides  criterion  for  accuracy. 

The  error  in  the  wave  velocity  for  long  waves  is  much  smaller 
for  quadratic  and  cubic  shape  functions,  even  if  one  takes  the 
fact  into  account  that  for  quadratic  and  cubic  shape  functions 
the  number  of  elements  needed  to  give  the  same  resolution  in 
Fourier  components  is  only  one-half  of  that  for  linear  wave 
functions.  (The  number  of  parameters  describing  the  solutions 
at  a  given  time  is  the  same  for  the  same  resolution.)  Short 
wave  errors  are  nearly  the  same  for  different  approaches.  In 
a  semidiscretized  method  one  is  led  to  a  system  of  ordinary 
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differential  equations  (in  implicit  form) .  The  solution 
must  proceed  in  small  steps  because  of  the  latent  presence 
of  particular  solutions  of  the  homogeneous  system  with  a 
short  wave  length.  Implicit  procedures  do  not  automatically 
insure  the  applicability  of  large  time  steps.  Short  waves 
(which  in  any  case  are  inaccurate)  can  be  suppressed  if  one 
works  with  third  degree  splines  rather  than  cubic  weight 
functions,  which  permit  discontinuities  in  the  second 
derivatives  at  the  grid  points.  Next  the  time  integration 
by  finite  elements  is  discussed.  (It  can  be  separated  from 
the  space  discretization  if  one  deals  with  bilinear,  biquadratic 
or  bicubic  elements.)  For  linear  elements  the  time  dependence 
can  be  treated  in  the  same  manner  as  the  space  dependence. 

For  Courant  number  1,  the  discretization  errors  for  time  and 
space  cancel  each  other.  The  method  becomes  unstable  if  the 
Courant  number  exceeds  one.  From  a  heuristic  point  of  view 
such  an  approach  is  somewhat  suspect,  for  it  does  not  quite 
fit  the  idea  of  a  marching  procedure  (which  can  always  be 
carried  out  in  hyperbolic  problems) .  An  analyogous  approach 
for  third  degree  elements  is  always  unstable.  A  modified 
approach  for  linear  elements  which  avoids  this  heuristic  dis¬ 
crepancy  gives  stable  but  strongly  damped  particular  solutions. 
From  the  point  of  accuracy,  one  would  therefore  need  rather 
small  elements.  A  similar  approach  for  third  degree  shape 
functions  converges  for  all  Courant  numbers  and  gives  tolerable 
results  up  to  Courant  number  one,  if  one  takes  constant  weight 
throughout  the  new  time  interval  in  combination  with  a  colloca¬ 
tion  condition  for  the  new  time  point.  An  example  discusses 
the  choice  of  the  mesh  if  ,  in  a  steady  flow  field,  the  Mach 
number  approaches  one  at  a  certain  line.  Even  then,  one  can 
always  find  a  mesh  which  guarantees  stability.  In  the  example 
under  discussion  the  particular  solutions  belonging  to  the 
highest  frequencies  (which  are  prone  to  cause  instability)  are 
of  importance  only  in  the  immediate  vicinity  of  the  sonic  line. 

A  rational  basis  for  the  choice  of  weight  functions  can  be 


obtained  by  relating  them  to  the  Green's  functions.  It  is 
possible  to  find  weight  functions  which  suppress  the  dominant 
long  distance  effect  of  a  residual.  For  the  Helmholtz  equa¬ 
tion  this  is  only  approximately  correct  and  requires  that  the 
elements  be  sufficiently  small.  For  the  hyperbolic  problem 
this  cannot  be  done  because  the  hyperbolic  distance  is  zero 
for  points  which  lie  on  the  same  characteristic.  A  better 
cancellation  of  long  range  effects  in  two-dimensional  problems 
can  be  obtained  if  one  chooses  characteristics  as  element 
boundaries.  The  weight  functions  should  then  be  constant 
in  one  of  the  characteristic  directions.  The  idea  cannot 
be  carried  over  directly  to  the  three-dimensional  problem. 
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SECTION  I 
INTRODUCTION 


This  report  discusses  the  application  of  the  finite 
element  concept  to  hyperbolic  problems.  Some  of  the 
observations  are  probably  familiar  to  those  who  have  been 
working  in  this  area;  nevertheless  the  author  hopes  that  the 
fairly  systematic  discussion  carried  out  here  has  some  value 
in  so  far  as  it  clearly  shows  certain  inherent  difficulties. 

An  attempt  to  develop  a  general  intuitive  background 
for  this  problem  area  will  be  made  in  Section  X.  In  the 
introduction  we  restrict  ourselves  to  one  general  observation 
which  motivates  the  choice  of  the  approaches  to  be  studied. 

In  the  classical  applications  of  the  finite  element  method  to 
elliptic  equations,  the  problem  is  governed  by  an  extremum 
principle.  This  gives  a  guarantee  of  convergence,  provided 
that  the  numerical  process  imitates  the  search  for  an 
extremum.  (Of  course  the  distance  definition  with  respect 
to  which  convergence  is  obtained  need  not  coincide  with  the 
error  characterization  desirable  for  an  engineering  application 
of  the  results.)  For  the  hyperbolic  problem  no  extremum 
formulation  exists.  In  the  classical  examples  one  can 
define  a  functional  which  is  stationary  if  oi>e  imposes 
suitable  boundary  conditions.  But  these  boundary  conditions 
are  not  identical  with  those  of  a  well  posed  hyperbolic 
problem.  For  practical  work  a  mere  variational  formulation 
has  only  limited  usefulness  to  begin  with  (because  of  the 
absence  of  guaranteed  convergence) ,  but  here  the  variational 
formulation  does  not  even  apply  to  the  problem  at  hand. 

Reddy  *  has  found  a  variational  formulation  for  some  hyperbolic 
equations  in  which  the  boundary  conditions  are  properly  taken 
into  account.  A  simple  example  is  shown  in  Appendix  I,  but 

*  (1)  Reddy,  J.  N.,  A  Note  on  Mixed  Variational  Principles 

for  Initial-Value  Problems,  Quarterly  Journal  Mech.  Appl. 

Math.,  Vol.  XXVIII,  Pt.  1,  1975,  pp’.  123-132. 
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because  of  the  lack  of  an  extremum  property  even  this  for¬ 
mulation  does  not  seem  useful  from  a  numerical  point  of  view. 

Sections  II  through  V  study,  for  the  wave  equation, 
different  approaches  by  semi-discretization;  that  is,  the 
finite  element  concept  is  applied  only  to  the  spatial 
direction.  This  leads  to  systems  of  ordinary  differential 
equations  with  the  time  as  independent  variable.  It  is 
assumed  that  this  system  is  solved  by  an  accurate  integration 
method.  In  Section  VI  the  finite  element  concept  is  applied 
also  in  the  time  direction.  This  discussion  becomes  fairly 
simple  because  of  the  special  choice  of  bilinear,  biquadratic 
or  bicubic  shape  functions.  The  results  show  the  importance 
of  the  Courant  number.  Sometimes  the  occurrence  of  large 
Courant  numbers  is  unavoidable.  Section  VII  gives  details 
for  a  simple  example  of  this  kind.  Section  VIII  tries  to 
provide  a  rationale  for  the  choice  of  the  weight  functions  by 
relating  these  to  the  Green's  function.  In  Section  IX  the 
finite  element  concept  is  combined  with  the  idea  of 
characteristics.  The  concluding  Section  X  develops  an 
intuitive  picture  of  the  problem  area.  In  addition,  it 
summarizes  the  specific  results  obtained. 

A  number  of  detailed  studies  are  included  in  the  form 
of  appendices.  We  already  mentioned  that  Appendix  I  gives 
an  example  for  the  variational  formulation  of  Reddy. 

Appendix  II,  following  an  idea  of  Soliman,  shows  that  rather 
general  boundary  conditions  (except  for  Dirichlet  conditions) 
fit  smoothly  into  the  concept  of  weighted  residuals.  Appendix 

III  gives  a  number  of  results  which  can  serve  to  check  certain 
formulae  which  occur  in  the  main  body  of  the  report.  Appendix 

IV  shows  that  the  study  of  integration  methods  (including  those 
based  on  the  finite  element  concept)  for  systems  of  ordinary 
differential  equations  can  be  reduced  to  the  discussion  of 

a  single  scalar  equation.  Appendix  V  discusses  Green's 
functions  for  the  wave  equation  in  two  and  three  dimensions. 


The  present  study  therefore  applies  the  finite  element 
concept  in  the  form  of  a  method  of  weighted  residuals.  In 
a  variational  formulation  the  weight  functions  which  are 
applied  to  the  residuals  are  always  identical  with  the  shape 
functions  used  to  represent  the  solutions  (if  not  directly, 
then  at  least  implicitly) .  In  the  absence  of  an  extremum 
formulation  this  choice  loses  its  usefulness.  We  shall 
admit  weight  functions  of  greater  generality.  (The  term 
shape  functions  will  be  restricted  to  expressions  which  serve 
to  represent  the  solutions,  the  word  weight  function  is 
self-explanatory.)  In  classical  applications  the  weight 
functions  have  finite  support.  They  are  or  for 

partial  differential  equations  of  order  two  or  four,  respectively. 
This  holds  in  particular  at  the  edges  of  the  region  of  support. 
This  latter  continuity  requirement  will  no  longer  be  imposed. 

In  the  integration  by  parts  needed  in  order  to  treat 
derivatives  of  generalized  functions  correctly,  certain 
additional  terms  outside  of  the  integrals  at  the  element 
boundaries  are  encountered .  With  the  traditional  choice 
of  the  weight  functions  (and  of  the  shape  functions),  these 
terms  will  vanish.  With  the  present  more  general  choice  of 
weight  functions  these  terms  must  be  taken  into  account. 

In  the  author's  opinion, this  is  a  technicality  which  should 
not  preclude  the  use  of  such  general  weight  functions. 

By  their  very  nature,  hyperbolic  equations  lend  themselves 
to  a  marching  procedure.  The  solution  at  a  certain  time  is 
not  influenced  by  what  happens  at  a  later  time.  We  have  limited 
our  discussion  to  approaches  which  can  be  interpreted  in  this 
sense.  The  procedures  obtained  in  this  manner  need  not  be 
stable.  A  major  part  of  the  present  work  is  the  study  of  the 
stability  (  and  of  the  accuracy)  of  different  methods  in  a 
simple  sample  problem.  In  most  of  the  report  we  deal  with 
rectangular  elements  and  bilinear,  biquadratic  or  bicubic 
shape  functions. 


SECTION  II 

EQUATIONS  OBTAINED  BY  SEMIDISCRETIZATION 

In  the  examples  of  this  report  we  shall  study  the  wave 
equation 

<b  —  <b  —  o 

yy  tt 

(If  one  replaces  <ftfc  by  $  then  one  obtains  the  linearized 
equation  of  two  dimensional  supersonic  flow  at  a  free  stream 
Mach  number  /? . )  Practically  in  all  cases  one  restricts 
stability  discussions  to  equations  with  constant  coefficients 
on  the  basis  of  the  argument  that  in  a  restricted  region  and 
for  a  very  fine  mesh  the  coefficients  are  practically  constant. 

The  choice  of  such  a  simple  equation  is  therefore  not  more 
restrictive  then  the  usual  stability  discussions.  We  assume 
that  the  region  extends  in  the  y  direction  from  to  +<». 

In  the  t  direction  the  solution  is  to  be  determined  for  a 

finite  interval.  In  this  section  an  approach  by  semidiscretization 

is  investigated;  the  differentiation  with  respect  to  the  t 

direction  is  retained  while  one  uses  a  finite  element 

representation  for  the  y  dependence.  Such  an  approach  is 

frequently  used  for  problems  with  two  or  three  space 

dimensions. 

In  this  section  we  shall  write  down  the  resulting 
systems  of  ordinary  differential  equations  for  different 
choices  of  the  (solution)  shape  functions  and  of  the  weight 
functions,  and  also  for  a  semidiscretized  finite  difference 
approach.  The  following  problems  will  be  treated. 

1.  The  y  derivative  is  replaced  by  a  finite  difference 
approximation . 

2.  0  is  approximated  by  a  piecewise  linear  function, 

<)>  is  continuous 

a.  the  weight  functions  are  identical  with  the  shape 
functions 
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b.  the  weight  functions  are  constant  over  intervals 
of  a  length  of  the  gridsize  which  straddles  the 
points  where  the  first  derivatives  are 
discontinuous 

3.  <t>  is  approximated  by  piecewise  third  order 
polynomials,  <j>  and  are  continuous 

a.  the  weight  functions  are  identical  with  the 
shape  functions 

b.  and  c.  the  weight  functions  are  constant  over 

intervals  of  a  length  equal  to  half  of  the 
grid  size.  Cases  b  and  c  differ  by  the  position 
of  these  intervals. 

4.  4>  is  approximated  by  piecewise  second  degree 
polynomials,  <J>  is  continuous. 

a.  the  weight  functions  are  related  to  the  shape 
functions 

b.  the  weight  functions  are  constant  over  intervals 
of  a  length  equal  to  half  of  the  grid  size. 

5.  <}>  is  approximated  by  third  degree  polynomials, 

4),  4>y,  and  <j>  are  continuous  (spline  approximation) 

a.  the  weight  functions  are  related  to  the  shape 
functions 

b.  the  weight  functions  are  constant  over  an 
interval  equal  to  the  grid  size. 

6.  is  approximated  by  piecewise  second  degree 
polynomials,  <j>  and  <j>y  are  continuous 

a.  linear  weight  function 

b.  constant  weight  function 

7.  <ji  is  approximated  by  third  degree  splines,  and  a 
collocation  method  is  applied. 

As  usual  the  solutions  are  represented  by  a  combination 
of  shape  functions,  defined  for  the  individual  intervals 
(elemental  shape  functions) .  The  formulae  arising  by  this 
standard  procedure  are  listed  below  for  the  convenience  of 
a  reader  who  wants  to  check  the  details.  They  are  readily 
tested  by  the  requirements  that  they  give  exact  representations 
of  the  operators  involved,  for  polynomials  of  a  degree  equal 
to  that  of  the  shape  functions.  This  is  done  in  Appendix  HI. 


The  interval  in  the  y  direction  is  h.  The  value  of 
y  at  the  kth  grid  point  is  denoted  by  y*.  For  the  vicinity 
of  the  point  y^  (usually  for  the  interval  y^_1  <  y  <  y^+1  we 
set 

Ay  =  y  -  yk  (2) 


The  elemental  shape  functions  are  written  as  functions  of 
(Ay/h) .  A  prime  denotes  the  derivative  with  respect  to 
(Ay/h) .  Furthermore,  let 

4„(t)  - 
^(t>  -  h 

The  elemental  shape  functions  are  given  by  the  following 
expressions. 

Linear  elemental  shape  functions  (Figure  1) 

M  f*y/h)  -  /  -  Ay /A 

(4 

K  (&y/h  )  •  aj /A 

They  satisfy 

A,  (o)  *  /  J  A/,0)  -  0 

AiCo)  =  o  ;  %(/)  *  / 


Quadratic  elemental  shape  functions  (Figure  2) 
A /(ay/A)  *  /  -  3  (Ay/h)  +  ifiy/hf 
A/Jdy/h)  *  -  (ay /A)  +  ^Oyftf 

A/3  (Ay/h  )  *  4 >(*y/h)  -tfry/hf 

They  satisfy 

a; Co)  -  /,  y('fi)  s  °,  VM*  0 

VJO)  *  %.(•/*■)  *  O  j  N^O)  -  f 

t^CO)  *  o,  a$0A)  *  '  j  %0)  =  o 


Third  degree  elemental  shape  functions  (Figure  3) 


1 

/vj*ylh)  - 

/VJ  fay/h)  *  (A>flh) 


“  J  ty/h)  4  l(A(f/h) 
3(ay/h)x  ~2.(Aylhf 
-  Ltylhf  4  fry/hf 
“  (*y/h)x  4  (*y, Jb) 


(6) 


They  satisfy 
, 

/I (Jo)=o  J 

W°)”°  , 
A>y  -  0 


*  0, 
^(0)  -  0, 
*  / 
*  0, 


>*/7y-  0 
/ 

y£tfy *  ' 

AJM  -  0, 

/i^/y  -  0 


yo)  *  0 
^V/;  -  0 
%'o)  *  0 
A^'y>y  *  ; 


With  these  characterizations  one  readily  obtains  the  following 
identities  which  can  be  used  for  checking  purposes. 

Linear  elemental  shape  function 

/  a  ^  f  ^4- 

<*y/4  - 

Quadratic  elemental  shape  functions 

/  -  4  -*4 
•ay/*  -  ^  »  A5A 
fy/h )*.  A i  */$/<■ 

Third  degree  shape  functions 

/  *  ^  *  A£, 

*y/A  r  A£  *  45  * 
f-fly  fh)***  4  L 

(aylh)  J-  %  *  JA£, 
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The  solutions  are  characterized: 

in  the  linear  case  by  the  values  of  <|>k(t), 

in  the  quadratic  case  by  the  values  of  4>k(t)  and  <J>k  +  jy2 

in  the  third  degree  case  by  the  values  of  4>k(t)  and 

The  values  of  k  range  in  our  examples  from  -»  to  +®  . 

Notice  that  one  has  two  parameters  per  interval  for  the 

quadratic  as  well  as  the  third  degree  case. 


(t). 


Let  w(Ay/h)  be  one  of  the  weight  functions.  Then  one 
obtains  an  approximating  equation 

J  &  wMb)d(^/h)  -  *  (7) 

For  each  choice  of  k  and  each  weight  function  we  obtain 
one  condition.  Here  one  must  substitute  for  <j>  the  approximations 
listed  above.  The  weight  functions  w  have  finite  support 
(that  is  they  are  different  from  zero  only  over  a  finite 
interval)  and  the  integrations  are  only  needed  for  finite 
intervals.  For  linear  and  quadratic  shape  functions  the 
term  will  lead  to  delta  functions  at  the  grid  points, 
the  integrals  must  then  be  evaluated  in  the  sense  of 
generalized  functions. 


Case  1.  Finite  difference  approximation. 
One  obtains 


"^1  0 


J 


-oo  <£,<**> 


(8) 


Case  2.  Linear  shape  functions 

The  shape  function  belonging  to  the  point  yk  is  given  by 
Figure  4. 

A't'oyAJ  *  *>JlAy/h)+>) 


(It  is  understood  that  the  elemental  shape  functions  N^  are 
zero,  unless  their  argument  lies  between  zero  and  one.) 
Accordingly,  N^  (Ay/h)  is  different  from  zero  only  for 
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yk  <  y  <  yk+1,  is  different  from  zero  only  for  yk-1  <  y  <  yk< 
One  has  in  the  interval  yk_^  <  y  <  y^^ 

fa  t)  -  ty((*ylh)+  •)  *  fat)  ft  ((Ayfhht)  *  Nt  Utfhf)  (g) 

NJy/k) 

Case  2a 

The  weight  function  belonging  to  point  yk  is  given  by 
(Figure  4) 

ufayjh )  *  N,  (*y!h  )  *  4  ( Ouflti  *  *) 

or 

»(Ayjh)  =  t+  fcyjh)  y  yi  f  <y  <  q 

v(aylh)  -  I  -  Uy/h)  y  ^  <y  <yiH 

Otherwise  they  are  zero. 

One  obtains 

£l  M  ♦  r<4k.,  -*i  '*W  ■ 0  HD 

For  a  test  see  Appendix  III. 

Case  2b 

The  weight  function  belonging  to  point  yk  is  given  by 
(Figure  5) 

*(*y/h)*  /  j  <  fay/h)  <  i  (i2) 

One  obtains 

jpM'  (13) 

For  a  test  see  Appendix  III. 


Case  3  Third  degree  elements 

The  shape  functions  belonging  to  point  yk  are  (Figure  6) 

(&y/h)  ■+  ^  faffrl) 

and 

4$  (&tflh)  *  4/f  ((Atf/h)  *  / 


In  the  interval  yk-1 


<  y  <  yk+1  one  has 


4 *>> f  tjjt) -J  fhiM  H> 

4  j \lt)[H((^lh*l)  *  ty/iy/h)]*  4  qrv/t)} 

+faW&y/h)  *faft)N¥!iylh)  (14) 


Case  3a 

The  weight  functions  belonging  to  point  y^  are  (Figure  6) 


yfay/Aj  -  ty(y/h)  f  Mtfty/hJ+f) 

and 

#L ( ty/h)  4/3 (&y I V  +  %■  (ty/h) */) 


(15a) 


(15b) 


One  obtains  the  following  differential  equations  (two  for 
every  value  of  k) 


felfah  fa -*i  *fa>  -  4h  fa,  - fa} 

+hLl- Tfa-ysfa)*  £ffa,  -A',)}  ‘° 


(16a) 


&{£,  fa- fa  fafafa,fafafafan 

/  *~X{-Wfa-A-, >*{■?&-  75fa-Wfa)l}-o 


j-"<i  <*> 


(16b) 


Case  3b 


The  weight  functions  are  given  by  (Figure  7) 


k/(4y/*tJ*/  /  -/  <  4j/A  <*•/* 


-/  i  -i<  itM  <° 

i  /  j  o  <  I1*  <  * 


(17a) 


(17b) 


otherwise  they  are  zero.  These  weight  functions  are  obviously 
equivalent  to 


v=/y  0  <4y/h  <±  and  /  y  -/ 


One  obtains 


j  -*>  ^4  <*>  (18a) 


+  ~lik  -&)}=  o 

* 1% " m rA'> 

+  h  /-  "  ° 

Case  3c 

The  weight  functions  are  given  by  (Figure  8) 

wjAy/h)^i  .  -  I  J  if  <  ay//,  <  i/m 

»i.(*ylk)m/  i  //*  <  4y/h  <*/+ 

otherwise  they  are  zero. 

One  obtains 

£  t.l  A  ^  J  .t U  I  _  _ /l  /  /'  •  > 


(18b) 


feUA 

**’V-  f  c J 


-*  <4  <* 


(20a) 
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and 


h  It  fa  +An)-  7 

-  to  f  a  <  a> 


(20b) 


Case  4 .  Second  degree  elements 

In  the  interval  yk-1  <  Y  <  Yk+1 
representation  for  4> 


one  has  the  following 


+  A  J*)  (&/>>*»  +<f.(*>Kftj/hj 

»'z  *>  +  JL  ' 

Case  4a 


The  weight  functions  are  given  by  (Figure  9) 

r  /  +  Uy/h)  /  -/  <  *yM  <  c 

l  /_ 


fay/h) 


o  <ay/4  <  * 


(21a) 


*i  (cyjh )  «  /V  (ay lit)  -  *H&y/h)  -  kfty/h)*' 


(21b) 


Notice  that  the  weight  function  w^  is  a  linear  combination  of 
elemental  quadratic  shape  functions.  This  choice  has  been 
made  to  avoid  negative  weights. 


One  obtains 

Te-j  7  A  -  ZA  +  Aft}  *  0 

-  *>  <  4  < 


(22a) 


and 

Tv-lfsWk'  +4+/J  *  /r  A*i\  +  h  7’jfy  *Ah)  +  I"  0 

-a*  <  4  <*> 


(22b) 
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Case  4b 


The  weight  functions  are  given  by  (Figure  8) 


•  )  b 

<  <  *  i 

(23a) 

,  i 

<  *9 ft  <  i 

(23b) 

One  obtains 

£■{£  A-  it  - 
-*(ah 

+  j-o°  <-£<*> 

(24a) 

(24b) 

-  to  <  4  <-  <to 

Case  5.  Third  degree  shape  functions  with  continuous  second 
derivatives 


The  jump  of  the  second  derivative  with  respect  to  Ay/h 
at  a  point  y  =  yk  is  given  by  the  left  side  of  the  following 
expression 


It  is  not  difficult  to  devise  a  test  for  the  correctness  of 
this  expression  similar  to  the  tests  derived  in  Appendix  III. 
This  expression  is  now  combined  with  Eqs.  (16a)  pertaining 
to  the  weight  function  eq.  (15a)  and  with  Eq.  (18a)  pertaining 
to  weight  functions  Eq.  (17a) .  Now  one  has  the  systems  of 
equations . 

Case  5a 

£  M-  *  h  fa-*  -  4 fa* '  A-/# 
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and 


6  tti*/  '  ~Z  fan 

-  oo  <  k,  <  a& 

Case  5b 

+  *‘Y-  r  (h-r*h  +fa") *  T  (h+*  -A-s}l 
£(fk+i  -A-*)  ~  *A' "  *{&,+&«) m  0 

-  00  <  k  <  &> 

The  weight  functions  of  the  case  3c  are  unsuited  to  such  an 
approach. 

In  practical  applications  there  are  always  a  finite 
number  of  values  and  <J>k'.  Then  one  can  express  the  values 
of  4»k  ’  in  terms  of  the  values  of  <p^  by  means  of  the  second 
equation.  One  thus  obtains  a  system  of  equations  for  the 
<J>k's  only.  For  the  present  discussions  this  is  not  practical. 

Case  6. 

Quadratic  shape  functions  studied  in  case  4  allow  the 
first  derivatives  to  be  discontinous .  It  might  be  a 
simplification  if  one  imposes  the  condition  of  continuity 
of  the  first  derivatives.  The  jump  of  the  first  derivative 
at  the  point  y  =  y^  is  given  by  the  left  hand  side  of  the 
following  equation 

~  6  A  “(h+l+A->)+  0  (28) 

Case  6a.  (linear  weight  functions) 

We  use  this  condition  in  combination  with  Eq.  (22a) 
which  holds  for  the  weight  functions  Eq.  (21a) ,  Figure  9a. 


(26b) 


(27a) 


(27b) 
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One  then  obtains 


*  *£ -i  ' h  f ~ZA  *  A*/ } 
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(29a) 


'  6  A  '  (h-*  y  A-h)  +  ¥Wk,+i 

-  OO  <  k  <  oo 


(29b) 


Case  6b.  Constant  weight  function 

The  expression  Eq.  (23a)  is  not  suitable  as  sole  weight 
function  because  it  fails  to  cover  the  whole  interval.  We  use, 
instead,  the  weight  Eq.  (17a)  (Fig.  7a)  and  obtain 

"  IvVoi'f+i-,)  (30a) 

*  fW  * 

-  id  -  rd  *4  )  ++M  , +4  ,)»  o  (30b) 


Case  7 . 

The  second  derivative  is  defined  everywhere  if  one 
deals  with  third  degree  polynomials  and  continuous  first 
and  second  derivatives.  It  is  then  possible  to  use  such 
an  approximation  in  order  to  define  the  second  derivatives 
at  the* grid  points  and  to  apply  a  collocation  method.  For 
4>  at  point  yk  one  obtains  from  the  interval  yk<  y  <  y^+^ 

y  (^ir/ "  J 

correspondingly 

" h 6 A  ¥  } 

Eq.  (28)  guarantees  that  these  expressions  be  the  same. 

For  the  sake  of  symmetry  we  write  ^.,v(y,,)  as  the  average 

yy  k 

of  these  values.  Then  one  obtains  from  Eq.  (1)  (with  a 
change  of  sign) 
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-  h  {sffk-t  B  0 

(31a) 


together  with 

6  ' &-•)  ‘  '  2  <&//  yii4  ^  -  ° 

-op  <  k,  <  oo  (31b) 
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SECTION  III 

STABILITY  TESTS  FOR  THE  SEMIDISCRETIZED  PROBLEM 


A  stability  analysis  can  be  carried  out  in  the  usual 
manner.  The  interval  in  the  y  direction  is  assumed  to  go 
from  -»  to  +».  Alternatively,  one  might  assume  periodicity 
conditions.  Then  one  sets 

^(6)  =  C  (ifskh)  *y(cyt)  (32a) 


and  if  needed 

«•  iC  (if'M)  *%(>  (ivt  )  (32b) 

and  computes  the  values  of  v  in  dependence  upon  y.  The  method 
is  stable,  if  v  is  real  or  has  a  positive  imaginary  part. 

If  y  is  small  (long  waves),  then  one  expects  that  v  =  +  y. 

In  a  precise  solution  of  the  original  problem  this  result 
should  be  found  for  all  values  of  y.  In  an  approximate 
method  this  relation  will  be  satisfied  only  approximately. 

This  gives  an  insight  into  the  errors  introduced  by 
different  approaches. 

One  obtains  with  Eqs.  (32) 

4-/  *  C  **/>  (•+)  (33) 


-  C  ty>(<:£»kh)M/3  (Cut) (iC 


&?/  -  *  C vy>(Cf,kh)  txp(iyt)f-Z)  ) 


(34) 

(35) 

(36) 
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One  obtains: 

Case  1.  (finite  differences,  Eq.  (8)) 


fh  £  tb  v  C^h/2. ) 


(36) 


Case  2a.  (linear  shape  functions,  linear  weight  functions 
Eq.  (11)) 

-  ***(!**&  (37) 

Case  2b.  (linear  shape  functions,  constant  shape  functions, 
Eq.  (13)) 

/-(t/zizM-fah/zl  l  ) 

Case  3a.  (third  degree  shape  functions,  weight  functions 
derived  from  shape  functions  Eqs  (16)) 


(z  ¥/s)  sc*i£uA/z )  -(ifc) 

(t/S )  (t/s)  +  f l/*S) Jivv'feAsh/i. ) 

l 

/  -  (li/is)  (te/zto)  ) 


(39) 


\  r 


(/3 /Z/O)  **  )  (t/z/O)  *  (//JS)A*%u4U) J 


-  O 


In  order  for  this  system  of  eauations  to  have  a  nontrivial 

solution  for  the  vector  [c,  cJ  r  its  determinant  must 

.  .  2  2 
vanish.  This  gives  a  quadratic  equation  for  v  h  .  For 

the  numerical  evaluation  this  determinant  is  a  suitable 

starting  point. 
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Case  3b.  (third  degree  shape  functions,  constant  weight 
functions  Eqs.  (18)) 


6  )  -  C/zJ  h 

-  i  jC*£*h)  3  ~ 


f 


-  v^h 


a- 


P//b)+ 

* 


(40) 
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Case  3c.  (third  degree  shape  functions,  constant  weight 
functions  Eqs.  (20)) 


(  ($/i  )  *i*t(f,h/z  )  | 

0 


i .z* 


f  (t/z)  -  {7//Z4)  firffcJi/l)  (/3//SiC)ji*£~i»)  1 


[  (t/z-  )cn(f,h/z) 


(ft/m  )/rn^uWzJ 


(41) 


Case  4 .  Quadratic  shape  functions 
We  set 

fa  -  C 

fa* l  »  Zu{>{Lf,(k,H)h)  **/>(***) 


19 


Case  4a.  (weight  functions  derived  from  shape  functions, 
Eqs.  (21)) 


«  O 


tf  Jt»f(gch/x  )  o 

'/j 

i 

tt/i 

(i/tf)eos(gJi/L) 

*/iC 

J 

(43) 


Case  4b.  (constant  weight  functions,  Eqs.  (23)) 


b  -4 '  ces^Jt/x) 

v'-h1 

(7/lif)  f-  {i//2.)fcHj^h/X  )  ( S/ltt)COi(£*.h/i) 

-  V  cos^h/x)  b 

-  v  n 

i/i  b  CM^uh/x)  n/x  b 

L 

y\ 

(44) 

Case  5.  (third  degree  shape  functions  with  continuous 
second  derivatives) 

Eq.  (25)  which  introduces  continuity  of  the  second  derivatives, 
combined  with  Eqs.  (32)  gives 


hi 


-  v/frsh/l ) 


(45) 


Case  5a.  (weight  function  related  to  shape  function) 

One  obtains  from  the  first  of  Eqs.  (39) 

0  -  Ot/iS)  h/x ))  (/- ))  ■*  fa/i/t )  &*Sfc**/ 

This  can  be  simplified  to 

1  ^  /i.  J  **  ^ij/r _ _ 

(/  -  *  fio/iOS)  )  (46) 

Case  5b.  (weight  function  constant) 

One  obtains  from  the  first  row  in  the  determinant  Eq.  (40) 


^  -  fa)  Jaffa**) 

(/  -  faft)  Jt+Ljguk/jj)//-  ■+  ) 


This  simplifies  to 

,  * icyC-(^)  -  2.Ji>J*(fJ*U)  (4?) 

’  /-  te/6]Ji%fy*A)  +  f'frj  > 

Case  6 .  (quadratic  shape  functions  with  continuous  first 
derivatives 

The  condition (28)  for  continuity  of  the  first  derivatives 
in  combination  with  Eq.  (42)  gives 

C'm.  - ~ P/*) £ 
toi  *>/*•) 

Case  6a.  (linear  weight  function) 

One  obtains  from  the  first  row  in  the  determinant  Eq.  (43) 
t/i.  b  ti+i? (/Ji/t ) 

»  *  »  - - -  (48) 

/  - 

Case  6b.  (constant  weight  function) 

The  result  is  derived  from  Eq.  (30a) 


b  ) 


(49) 


Case  7 .  (third  degree  splines  combined  with  collocation) 

Here  the  condition  (45)  (continuity  of  second  derivatives) 
applies  again.  One  obtains  from  Eq.  Ola) 

/  -  (i/d  )  / - /i/jJJtHtyY*}  (50) 
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SECTION  IV 

DISCUSSION  OF  THE  RELATIONS  DERIVED  IN  SECTION  III 

The  relations  between  vh  and  yh  derived  in  the  preceding 
section  have  been  evaluated  numerically.  The  information 
so  obtained  is  valid  only  if  one  solves  the  systems  of 
ordinary  differential  equations  which  are  written  down  in 
Section  II  precisely.  In  reality,  one  must  apply  some 
approximation  also  for  this  purpose.  Some  observations  regarding 
such  integrations  are  made  in  the  next  F«*ction.  The  formulae 
in  Section  III  are  based  on  the  hypothesis  (32a)  (supplemented 
by  other  expressions  of  a  similar  form)  that  is 

*  C  IX/3  (if*  k,fi  )■**/>(  1  1 

which  corresponds  to  an  approximate  particular  solution 

-  C  **/>  **/>( ±l»t)  -  ±  (vfcjt) 

The  wave  speed  in  this  approximation  is  then  +  v/y,  the  wave 
speed  for  exact  solutions  of  the  original  differential 
equation  is  +  1. 

In  the  approximations  discussed  so  far,  the  wave  speeds 
are  real  for  all  values  of  y;  there  is  no  damping.  (One 
would  speak  of  a  dispersive  method  because  waves  with  different 
wave  lengths  do  not  remain  together.)  Because  of  the  error 
in  the  wave  speeds,  one  obtains  a  phase  error  in  the  wave 
after  some  time  has  elapsed,  which  is  given  by 

(/  “  (' i 

In  the  graphs  the  value  of  v/y  is  shown,  the  phase  error 
decreases  more  rapidly  with  y  than  v/y  -  1,  because  of  the 
factor  y  which  occurs  in  the  last  expression.  The  contribution 
of  a  certain  wave  length  to  the  solution  becomes  meaningless 
if  the  phase  error  exceeds  some  fairly  small  number  (perhaps 
ir/6  or  smaller)  .  One  must  make  the  grid  size  h  small 


i 


enough  so  that  for  the  time  for  which  the  solution  is 
required,  the  phase  error  for  the  significant  wave  lengths 
stays  small. 

For  the  cases  1  and  2,  one  obtains  only  one  value  of 

2 

(vh)  for  each  choice  of  (yh)  .  For  cases  3  and  4,  one 

2 

obtains  a  quadratic  equation  for  (vh)  .  This  requires  an 
explanation. 

Consider  expressions  which  assume j at  the  grid  points^ 
the  following  values: 

or  A  -  (5 

^  a  C  ) 

These  expressions  can  be  rewritten  in  the  form 
^  «  £  CO\f(^h  "  C  u>\  +  UP*)  ) 

(j;^  -  ~  C  +***')&) 


Except  for  a  trivial  change  of  sign  (in  the  sine  function) , 
one  obtains  the  same  values  at  the  grid  points  y  =  kh,  if 
one  replaces  y  by  -y  and  or  changes  yh  by  a  multiple  of 
2ir.  The  above  expressions  can  originate  from  waves  of 
the  fotm 


C  vj/fry) 

« •*  c  y  ^ 


where  y^  *  y  mod  y  ,  or  y^  s  -y  mod 

This  is  the  well  known  phenomenon  of  aliasing.  For  the 
evaluation  of  the  formulae  of  the  previous  section,  it  is 
therefore  sufficient  if  0  <  y  <  t. 
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In  general  one  will,  of  course,  associate  the  values 
of  vh  obtained  from  these  formulae  with  the  smallest  possible 
values  of  y  which  are  admissible  according  to  the  above 
formulae.  But  in  principle,  there  is  no  reason  to  disregard 
other  values  of  y^  right  from  the  start.  Of  course  one  will 
consider  only  values  of  y  for  which  the  pertinent  values 
of  v  give  tolerable  approximations  to  the  wave  speed. 


In  evaluating  vh  for  quadratic  and  third  order  shape 

2  2 

functions,  one  obtains  two  values  of  v  h  for  one  value  of 
yh.  One  will  surmise  that  one  of  these  values  vh  belongs 
to  yh,  while  the  second  one  belongs  to  2ir-yh.  This  is  borne 
out  by  the  graphs.  For  an  analytical  explanation  consider 
the  case  y  =  0.  Then  one  has  for  y^  =  y  =  0 


h  *"&)>  * 


/ 


hence 

and  for  **  (i-K/hJ 

<f>  - 


hence 

i £  ~  Juns  ((lZ/A)As/i )  m  O 

*  It  cat  ((**/Ukk)  -  / 


The  first  expression  ^  *  1  is  exactly  represented  by 
*  1  which  is  a  third  degree  polynomial.  The  second 
expression  gives  within  one  interval  the  expression 

j  s  It  ^  y  y¥  )  s  zc  (ay/h )((  -  0  -  ifayft)) 

It  is  shown  in  Figure  10,  together  with  the  function 
sin  (2iry)  . 
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This  result  is  not  surprising.  For  second  and  third 
order  degree  shape  functions,  one  has  twice  as  many  parameters 
per  grid  point  than  for  first  degree  shape  functions. 
Accordingly,  one  can  approximate  functions  with  a  wave 
length  one  half  of  that  which  can  be  approximated  by  first 
degree  shape  functions. 

The  amount  of  numerical  work  to  solve  a  system  of 
equations  with  double  the  number  of  unknowns  is,  of  course, 
at  least  twice  as  large.  To  compensate  for  this  fact,  one 
will  take  the  interval  in  the  y  direction  for  second  and 
third  order  shape  functions  twice  as  large  as  for  linear 
shape  functions.  The  scale  of  yh  for  second  degree  shape 
functions  has  therefore  been  chosen  1/2  of  the  yh  scale  for 
first  degree  shape  functions. 

Figure  11  gives  the  results  for  v/y  versus  yh  for  the 
cases  1  and  2.  Curve  a  pertains  to  a  finite  difference 
semidiscretization,  curve  b  belongs  to  case  2a  (linear  shape 
functions,  linear  weight  functions) ,  curve  c  to  case  2b 
(linear  weight  functions,  constant  shape  functions) .  Ideally, 
one  should  obtain  1  for  all  values  of  yh.  The  finite  element 
approximations  are  somewhat  better  than  the  finite  difference 
approximation,  the  region  of  values  yh  for  which  the  v/y  is 
nearly  correct  is  rather  limited.  Incidentally,  curve  b 
is  also  obtained  in  case  7  (third  order  shape  functions  with 
a  collocation  method) . 

Figures  12  and  13  give  results  for  third  degree  shape 
functions.  Notice  that  in  these  cases  yh  ranges  from  zero 
to  2tt  and  that  at  yh  =  tt  one  has  a  break.  The  reason  is 
explained  above.  Curve  12a  and  13a  give  the  result  for 
case  3a,  weight  functions  derived  from  shape  functions,  curve 
12b  shows  the  result  for  case  3b,  curve  13b  gives  the  result 
for  case  3c,  both  have  constant  weights  but  they  are  shifted 
with  respect  to  the  grid.  In  this  case,  weight  functions 


25 


1 


derived  from  shape  functions  are  superior.  It  is  rather 
disconcerting  that  a  small  change  of  the  weight  function 
as  it  occurs  between  Figures  12b  and  13b,  gives  a  rather 
significant  difference  in  the  results. 

It  is  conceivable  that  under  slightly  changed 
circumstances  the  results  would  be  different.  Under  the 
present  circumstances  the  procedure  of  case  3a  gives  the 
best  results,  but  the  author  is  not  sure  whether  or  not  the 
same  behavior  can  be  expected  for  more  complicated  problems. 

We  explained  above,  why  the  region  0  <  yh  <  2tt  for  the  third 
degree  shape  functions  is  considered  as  equivalent  to  the 
region  0  <  yh  <  tt  for  linear  functions.  For  shorter  wave 
lengths  the  error  in  the  wave  speed  is  about  the  same  for 
linear  and  third  order  shape  functions.  However,  the  values 
of  yh  for  which  long  waves  are  adequately  represented,  that 
is  where  v/y  is  close  to  1  is  considerably  enlarged  for 
third  degree  elements. 

Summarising,  one  can  say  by  taking  third  degree  shape 
functions  with  twice  the  grid  size  instead  of  first  order 
shape  functions,  one  is  able  to  represent  in  both  cases  the 
same  waviness  with  respect  to  the  y  direction.  The  wave 
speed  for  short  waves  is  falsified  by  about  the  same  amount 
in  either  case.  Long  waves  are  less  falsified  for  third 
degree  elements  than  for  first  degree  elements. 

If  one  uses  such  a  procedure,  one  must  always  take 
the  mesh  fine  enough,  so  that  the  behavior  of  all  the  waves 
that  are  important  for  the  solution  is  adequately  represented. 
Third  degree  order  elements  then  make  it  possible  to  use  a 
far  larger  mesh  (ultimately  to  work  with  fewer  unknowns) 
than  first  degree  elements.  Because  of  rounding  errors,  one 
must  expect  shorter  wave  lengths  to  appear  in  the  solutions,  the 
propagation  of  short  waves  is  falsified  by  the  procedure.  Such 
errors  do  not  die  out  (as  they  would  in  an  elliptic  problem)  but 
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at  least  their  amplitude  does  not  increase.  Moreover, 
if  such  waves  are  due  to  rounding  errors,  then  they  are 
likely  to  have  a  random  character  so  that  they  cancel  in  the 
average.  The  discretization  error  is  already  taken  into 
account  in  the  falsification  of  the  wave  speed. 

Since  short  waves  are  wrongly  represented  for  third 
degree  shape  functions,  it  might  be  desirable  to  suppress 
them  altogether.  This  can  be  done  by  replacing  one  set 
of  conditions  which  arise  by  the  use  of  a  certain  weight 
function  by  the  requirement  that  the  second  derivative  be 
continuous.  In  this  manner  one  connects  the  ’  s  with  the 
4)^'s.  In  other  words,  one  uses  a  third  order  spline 
representation.  The  results  are  shown  in  Figures  14  and  15. 

The  requirement  of  continuity  of  the  second  derivative  is 
a  restriction  of  the  space  of  functions  that  is  available 
for  the  representation  of  the  solution,  therefore,  one 
expects  some  deterioration  in  the  wave  speed.  This  is  indeed 
the  case.  Figure  14  compares  the  case  3a  (third  degree 
shape  functions,  weight  functions  derived  from  the  shape 
function)  with  a  case  5a  where  one  of  the  weight  functions 
is  replaced  by  the  requirement  of  continuity  of  the  second 
derivative.  Figure  15  gives  the  same  comparison  for  the 
case  3b  constant  weight  functions  with  case  5b.  One  has 
indeed  a  deterioration  of  the  result  for  shorter  waves. 

(In  these  curves  the  maximum  value  of  yh  is  it.  A  comparison 
with  Figure  11  shows  that  this  procedure  is  indeed  much 
better  than  for  first  degree  shape  functions. 

In  third  degree  splines,  the  second  derivative  is  defined 
everywhere.  It  is  then  possible  to  use  a  collocation 
method,  in  which  the  differential  equation  is  satisfied  for 
the  grid  points  y  =  kh  (case  7) .  The  results  agree  with 
curve  c  in  Figure  1.  A  comparison  of  different  third 
degree  formulae  is  shown  in  Figure  16.  Cases  5a  and  5b 
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(third  degree  splines  give  nearly  identical  results  which 
are  fairly  acceptable)  .  A  collocation  method  based  on  third 
degree  splines  is  not  an  improvement  over  first  degree 
functions . 

The  results  for  second  degree  shape  functions,  cases 
4a  and  4b,  are  shown  in  Figure  17.  In  Figures  18  and  19 
the  cases  4a  and  4b  (second  degree  shape  functions)  are 
compared  with  cases  6a  and  6b  (second  degree  splines) . 

Second  degree  splines  are  very  inferior. 

Figure  20  shows  a  comparison  between  the  results  for 
third  degree  shape  functions  (case  3b)  and  a  quadratic  shape 
function  (case  6b) .  The  comparison  is  between  a  moderately 
good  case  for  third  degree  shape  functions  and  the  best 
case  for  second  degree  shape  functions.  In  the  overall 
picture  the  two  methods  are  about  equivalent,  for  short 
waves  the  second  degree  method  gives  better  results.  However, 
for  long  waves  which  are  the  more  important  ones,  the  third 
degree  method  gives  a  better  approximation  to  the  wave  speed 
(v/u  =  1)  over  a  wider  range  of  values  ph.  The  amount  of 
labor  in  solving  the  system  of  ordinary  differential  equations 
is  the  same  for  third  and  second  degree  shape  functions. 

One  concludes  that  the  use  of  third  degree  shape  functions, 
or  the  use  of  third  degree  splines,  is  preferable.  However, 
this  conclusion  is  not  generally  correct  because  it  disregards 
the  labor  to  set  up  the  system.  This  seems  to  hold  in  particular 
for  elliptic  equations. 


SECTION  V 


REMARKS  ABOUT  THE  NUMERICAL  SOLUTION  OF  THE  SYSTEMS  OF 
DIFFERENTIAL  EQUATIONS  DERIVED  IN  SECTION  II 

The  discussion  of  Section  IV  presuppose  that  one 
solves  the  systems  of  differential  equations  which  arise  by 
various  methods  of  semidiscretization  perfectly.  In  reality, 
one  applies  also  for  this  purpose  some  numerical  approach. 

We  shall  see  that  the  error  incurred  in  the  solution  of  the 
differential  equations  for  the  time  dependence  combine 
with  the  errors  due  to  the  semidiscretization.  Under 
favorable  conditions  they  may  cancel  each  other.  These 
phenomena  will  be  discussed  later.  At  the  moment  we  discuss 
the  application  of  finite  difference  techniques  to  the  systems 
of  differential  equation  s  of  Section  II.  One  might  think  of 
some  carefully  written  code  with  provisions  for  error  control. 
Most  of  these  methods  make  use  of  information  generated  at 
several  preceeding  points  in  time. 

One  notices  that  all  formulae  based  on  the  finite 

element  approach  fail  to  give  explicit  expressions  for  the 

2  2  ?.  2 

time  derivatives  d  <J>k/dt  or  d  ’0^/dt  .  To  obtain  these 
derivatives  in  terms  of  the  value  of  the  functions  <J>^  and 
()>£,  one  must  solve  a  system  of  linear  equations.  For  one 
space  dimension  this  task  is  not  too  time  consuming.  One 
then  obtains  a  banded  matrix  and  the  system  of  equations  can 
be  solved  very  efficiently.  The  problem  becomes  cumbersome  for 
more  than  one  space  dimension.  For  two  space  dimensions,  one 
obtains  a  block  tridiagonal  matrix  where  the  size  of  the 
blocks  depends  upon  the  number  of  grid  points  in  one  space 
direction,  thus  the  individual  matrices  may  be  rather  large. 

An  iterative  solution  of  these  systems  may  be  possible  in 
the  cases  1  and  2  (Eqs.  (11)  and  (13))  and  also  in  Eqs.  (37) 
and  (38).  This  amounts  to  replacing  the  denominators  (which 
come  from  the  time  derivatives)  by  the  series 
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The  series  converges  for  the  values  of  a  which  occur  in 
these  formulae.  The  author  has  not  explored  the  question 
for  more  than  one  space  dimension. 

One  is  tempted  to  use  implicit  methods  for  the  solution 
of  the  differential  equations.  Implicit  methods  are 
particularly  useful  for  stiff  differential  equations;  that 
is,  for  equations  in  which  the  matrix  governing  the 
linearized  system  (after  it  has  been  resolved  with  respect 
to  the  derivatives)  has  eigenvalues  with  large  negative  real 
parts.  If  one  solves  such  systems  by  the  Runge  Kutta  method, 
or  one  of  the  predictor  corrector  methods,  then  the  step  in 
the  t  direction  is  limited  for  reasons  of  stability  to  about 
the  reciprocal  of  the  largest  eigenvalue.  A  suitable  implicit 
method  removes  this  limitation.  The  step  which  is  automatically 
chosen  by  the  routine  on  the  basis  of  an  accuracy  check 
is  determined  by  the  requirement  that  within  the  current 
step  the  solution  can  be  represented  with  a  desired  accuracy 
by  a  polynomial  of  a  chosen  order.  This  allows  the  method  to 
proceed  in  large  steps  in  regions  where  the  solution  is 
smooth.  Unfortunately,  these  preconditions  are  not  met  under 
the  present  circumstances.  The  eigenvalues  are  imaginary. 

One  cannot  expect  to  obtain  smooth  functions  <j>^(t)  since 
the  eigensolutions  pertaining  to  large  eigenvalues  (short 
waves)  do  not  die  out. 

Here  one  might  make  the  following  argument:  If  a  small 

step  is  needed  in  the  y  direction  to  express  the  solution 

with  sufficient  accuracy,  then  one  needs  in  principle  a 

corresponding  small  step  in  the  t  direction.  Otherwise,  one 

will  disregard  information  which  has  been  considered  as 
i 
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essential  in  the  initial  conditions.  From  the  point  of  view 
of  accuracy,  one  cannot  permit  in  implicit  methods  a  time 
step  which  is  much  larger  than  the  one  permissible  in  explicit 
procedures,  but  the  same  limitation  will  probably  be 
encountered  also  from  the  point  of  view  of  stability  of 
the  numerical  integration  procedure. 

The  implicit  method  requires  the  solution  of  a  linear 
system  of  equations  at  each  time  step,  even  if  the  differential 
equation  gives  the  derivatives  explicitly.  But,  under  these 
circumstances  it  is  not  necessary  to  have  the  derivatives 
in  an  explicit  form.  To  show  this  in  a  schematic  manner, 
let  us  consider  the  system  of  first  order  equations 

/y  ^  My  =  r 

where  y  and  r  are  vectors  and  L  and  M  are  matrices.  In  an 
implicit  scheme  one  ultimately  arrives  at  an  equation 

yfti-h  )  +  h  yYt-tfiJm  u 

where  u  is  a  known  vector  determined  from  the  values  of  y 
and  y* at  preceeding  points  in  time  and  on  the  value  of  r(t). 

The  specific  form  of  the  equation  depends  upon  the  integration 
method.  One  obtains  by  multiplying  this  equation  from  the 
left  with  L 

Ly(t+h  )  4  arwt  h  Lfft+t )  •  Lit 


(L  -  cinsih  M) )  *  Ik  -  censi  h  r* 

In  order  to  obtain  y(t  +  h) ,  one  must  therefore  solve 
a  system  with  the  matrix  L  -  const,  h  M.  If  the  derivatives 
are  explicitly  available,  then  L  is  replaced  by  the  identity 
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matrix.  If  one  wants  to  use  this  idea,  it  may  be  necessary 
to  rewrite  existing  routines  for  the  implicit  integration  of 
differential  equations. 

According  to  an  idea  of  Soliman  et  al.,  the  inversion 
of  the  operator  L  in  multidimensional  problems  is  facilitated 
if  the  weight  functions  and  the  shape  functions  have  a 
special  form.  In  the  two  dimensional  problems,  one  will 
introduce  a  two-dimensional  grid  (say  a  rectangular  grid) 
with  grid  points  characterized  by  two  subscripts.  To  the 
point  i,k  there  belongs  a  shape  function  N...  and  also  a  weight 

X ,  K 

function.  Assume  that  the  shape  function  and  the  weight 
functions  appear  in  the  form 

The  inversion  of  the  matrix  L  can  then  be  reduced  to  the 
repeated  solutions  of  one-dimensional  problems.  The 
price  for  this  simplification  is  the  restriction 
in  the  choice  of  shape  and  weight  functions.  Practically 
only  very  simple  shapes  are  possible;  for  instance, 
piecewise  linear  functions  in  a  rectangular  grid.  Even 
then,  a  direction  inversion  of  the  matrix  (L  -  const  hM)  is 
not  possible.  If  the  second  term  is  small  enough,  an 
additional  iteration  would  be  applicable.  This  restriction 
of  flexibility  (which  to  some  extent  contradicts  the  basic 
philosophy  of  the  finite  element  approach)  may  be  worthwhile 
for  multidimensional  problems  where  the  solution  of  a  large 
system  of  equations  is  a  very  time  consuming  element. 


SECTION  VI 
TIME  DEPENDENCE 

In  the  preceding  sections  the  original  partial  differ¬ 
ential  equation  has  been  reduced  to  a  system  of  ordinary 
differential  equations  (in  our  particular  example  with 
constant  coefficients) .  Appendix  IV  derives  the  familiar 
fact  that  the  results  of  the  treatment  of  this  system  by 
finite  difference  or  finite  element  methods  remains  unchanged 
if  one  first  makes  a  transformation  which  brings  the  system 
into  its  diagonal  form.  It  is  therefore  possible  to  con¬ 
sider  one  component  at  a  time.  The  following  discussions 
are  therefore  restricted  to  the  equation 

+VxU,s  o  (54) 

where  u  is  a  scalar  quantity;  the  constant  v  stands  for  one 
of  the  values  of  v  computed  in  Section  III  as  a  function  of 
the  reciprocal  of  the  wave  lenght  (except  for  a  factor  of  2tt)  . 
In  practice  one  does  not  make  such  a  diagonalizing  transforma¬ 
tion,  the  values  of  v  never  appear,  and  the  shape  functions 
and  weight  functions  used  in  solving  the  system  are  independent 
of  v. 

In  Section  II  shape  functions  for  the  y  direction  have 
been  used.  The  parameters  upon  which  they  depend  are  functions 
of  time,  which  are  now  represented  individually  by  shape 
functions  of  the  same  character.  The  scalars  u  in  Eq.  (54) 
are  linear  combinations  of  such  parameters.  One  arrives  at 
shape  functions  in  the  y,t  plane  of  considerable  complexity 
in  spite  of  the  fact  that  the  discussion  is  carried  out  for 
one-dimensional  problems.  Familiar  forms  are  bilinear, 
biquadratic  and  bicubic  shape  and  weight  functions  in  a 
rectangular  grid. 
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We  distinguish  between  steps  h  and  h.  in  the 

S  t 

space  and  time  directions  respectively.  The  stability  of 

2  2  2 

the  time  integration  depends  upon  v  hfc  .  The  maximum  of  v 
which  occurs  in  a  specific  problem  depends  upon  the  step  h 

s 

in  the  space  direction.  Sometimes  it  may  be  desirable  to 
have  a  method  which  is  stable  even  if  vh^  is  large. 

Hyperbolic  problems,  by  their  very  nature,  are  initial 
value  problems?  at  an  initial  value  of  t  the  values  of  u  and 
du/dt  in  Eq.  (54)  are  prescribed.  The  solutions  can  there¬ 
fore  be  obtained  by  a  marching  procedure.  We  shall  restrict 
ourselves  to  methods  which  can  be  interpreted  in  this 
manner . 

The  application  of  finite  difference  methods  has  been 
discussed  in  general  terms  in  Section  IV.  Here  we  discuss 
a  specific  approach  which  is  sometimes  advocated.  One  writes 
Eq.  (54)  as  a  system  of  first  order  equations 


tju, 

di t 


-  Ur 


du*  + 

dt 


0 

o 


(55) 


We  use  dots  to  indicate  differentiation  with  respect  to  time. 

In  Section  II  differentiations  with  respect  to  y/h  have  been 

s 

denoted  by  prime. 

We  use  equidistant  points  in  time 


t 


k 


kh 


t 


and  set 
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Let  0  be  a  constant;  0  <  0  <  1.  Then  one  obtains  in  a  familiar 
manner 


- v -4  *(>-»*&’  ° 
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From  this  system  one  finds  uk+1  and  uk’+1  in  terms  of  and 
UjJ,  To  study  the  stability  one  introduces  an  amplification 
factor  p. 


(57) 
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This  leads  to  the  equation 

[9  + 


-ht[9+?O-0)] 

-/+f 


(58) 


Hence  (f  *  fO- 0 

/  i  C  0 
/  ft  *  (»*!£){'- 6) 
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One  has 


/,/  J  1 


<4 


(59) 


The  method  is  stable  for  6  «  1/2.  One  also  has 

tyf  •  i  *  ortffyfi-u)}  (60) 

A  single  wave,  defined  by  its  value  at  the  grid  points  is 
given  by 

f  -  COS0*y) 

»  l  ■**/*( \  ^  SyfJf 

WtOl  cot  fa  ^-[atciy(^)  t  arc$fi\O-6J0} 

1  (61) 

If  there  were  no  phase  error  the  factor  of  t^  in  the  last 
term,  namely 

-i  [<uctf  (Vi.*l  +  *rct}(r\('-e»] 
would  be  equal  to  y. 

Depending  upon  the  specific  form  of  the  differential 
equations  derived  in  Section  II  and  upon  the  choice  of  hg 
one  has  different  relations  between  v  and  y.  Ultimately, 
one  is  interested  in  the  relation  between  y  and  quantities 
characterizing  the  amplification  factor.  For  some  values 
of  h  /h.  and  for  some  values  of  y  a  cancellation  of  the 

S  L 

discretization  errors  may  occur. 

mhe  relations  between  y  and  v  for  cases  1  and  2,  Eqs. 
(36),  (37),  and  (38)  can  be  written  as: 
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with  0=0,  2/3,  and  1/2,  respectively.  For  these  values 
of  a  and  for  different  values  of  the  Courant  number  $  = 
ht/hg,  the  following  expression  is  shown  in  Figs.  21,  22,  and 
23. 


<r 


This  expression  arises  from  Eq.  (61)  for  3  =  1/2.  In  this 
case  I p |  =1  and  one  has  no  damping.  Ideally  this  expression 
ought  to  be  1  for  all  values  of  u  .  The  phase  errors  which 
arise  are  by  no  means  small.  The  curves  shown  in  Fig.  11 
supplement  these  figures.  They  give  the  results  for  Courant 
number  0.  In  all  cases  Courant  number  1/2  gives  about  the 
best  results,  particularly  in  conjunction  with  method  2b. 

The  results  deteriorate  considerably  with  increasing  Courant 
number. 


Next  we  study  cases  in  which  the  discretization  applied 
in  Section  II  for  the  space  direction  are  used  also  for  the 
time  direction.  We  restrict  ourselves  to  cases  1  through  3. 
In  actual  computations  4>  ( y^,  and  in  cases  3 

also  and  t^)  are  known.  Here  t^  refers  to 

the  time  for  which  the  computation  has  just  been  completed, 
l  ranges  over  all  stations  in  the  y-direction.  The  procedure 
then  computed  4*(y^»tjc+^)  and  in  cases  3,  also  <|>Ty^  ,tk+1) 
According  to  the  observations  made  at  the  beginning  of  this 
section  it  suffices  if  one  studies  Eq.  (54).  Then  one 


determines  u^+1  from  ^  and  uk  and,  in  cases  3,  u^+1  and 


uk+l  from  uk_1,  u^,  and  u^.  One  obtains  indeed  a 
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inarching  procedure. 

Case  2a  gives 

»"[%  *  £  (%rl%  *  4  i  (%h  ~iui.  4\J ' 0 

In  the  present  context  v  is  given  as  a  function  of  yhg. 

The  formulae  are  found  in  Section  III.  Particular  solutions 
of  this  equation  are  found  by  the  hypothesis 

tOA  m  JXfi  fifty) 

where  is  unknown.  One  then  obtains 

*)-o 


Hence,  for  case  2a,  after  substitution  of  Eq.  (37) 

.V  .i  ***(££)  _  » tin  (~> 

For  h./h  =  1  (Courant  number  1)  one  obviously  has 
t  s 

*J  -  if 

This  means  that  the  wave  speed  is  correct.  The  errors  due 
to  the  discretization  for  the  t  direction  cancel  the  errors 
of  the  y  discretization.  In  the  interval  0  <  v  h.  <tt  the 
expression  on  the  right  of  Eq.  (62)  is  a  monotonically 
increasing  function  of  sin  (v^ht/2) ,  its  maxium  is  12. 

If  the  left  hand  side  exceeds  this  value,  then  the  value  of 
for  which  this  equation  is  solved  cannot  be  real.  From 
the  two  conjugate  complex  solutions  which  will  then  appear, 
one  will  give  waves  whose  amplitudes  increase  with  time. 
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For  h,_/h  <  1  one  has  a  real  value  of  v,  for  each  choice 

t  s  1 

of  u;  the  method  is  stable  for  Courant  numbers  smaller  than 
1.  The  cancellation  of  discretization  errors  for  Courant 
number  1  has  no  counterpart  in  elliptic  equations. 

The  procedure  for  third  degree  shape  functions  can  be 
carried  out  in  analogous  manner.  However,  one  encounters  a 
significant  difference.  Consider  case  3a.  The  relation 
between  y  and  v  is  given  by  Eq.  (39) .  It  has  been  obtained 
from  Eq.  (16)  by  setting 

W  =  *k  exp  (iwhs) * 

For  the  present  discussion  is  it  preferable  if  we  write 


*k+i  "  V 


(63) 


where  p  is  some  point  on  the  unit  circle  in  the  complex  plane. 
Then  one  obtains  instead  of  Eq.  (39) 


-jrff-l  tf') 
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For  p  on  the  unit  circle,  this  equation  gives  real  values  of 
2 

v  h  .  Now  we  apply  a  discretization  analogous  to  case  3a  to 

S 

Eq.  (54) .  The  results  can  be  written  down,  by  making  appro¬ 
priate  modifications  in  Eq.  (16). 
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(67) 


These  values  of  p,  lie  on  the  unit  circle  in  the  complex 

2 

plane.  The  value  of  v  pertaining  to  these  solutions  is 
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But  there  exists  a  second  solution  for  v2 .  We  shall  see 
that  among  values  of  belonging  to  this  second  solution, 
there  is  always  one  which  gives  an  unstable  solution.  It 
follows  from  Eq.  (65)  that  the  product  of  the  two  roots  for 
P1  is  alwaYs  1 •  All  nonincreasing  particular  solutions 
must  have  values  of  pL  lying  on  the  unit  circle,  for  there 
belongs  to  each  root  which  lies  inside  the  unit  circle 
another  one  outside.  If  one  chooses  p1  as  a  point  of  the 
unit  circle  then  one  can  evaluate  Eq.  (65)  for  v2hfc2. 

(This  is  the  evaluation  which  led  to  Fig.  12  curve  a,  but 
now  we  consider  phj.  and  vhfc  as  coordinates.)  To  each 
value  of  yhfc  one  obtains  two  values  of  vht,  each  a  monotonic 
function  of  yh^;  these  two  functions  do  not  overlap,  it 
follows  that  for  these  curves  (which  exhaust  all  possibilities 
for  p^  being  on  the  unit  circle)  one  obtains  for  each  choice 
of  vhfc  one  value  of  yhg  (and  consequently  one  value  of  v2) . 

The  slope  of  these  curves  does  not  vanish,  accordingly  there 
are  no  double  roots.  It  follows  that  the  second  root  for  v2 
gives  values  of  p^  off  the  unit  circle,  one  of  them  pertains 
to  an  unstable  solution.  The  same  argument  can,  of  course, 
be  made  if  h./h  ^  1. 

t  S 

This  observation  can  be  connected  to  a  phenomenon  which 
arises  if  one  solves  ordinary  differential  equations  by 
integration  techniques  in  which  one  uses  information  generated 
for  preceding  points  in  time.  One  then  deals  in  principle 


with  difference  operators  which  are  of  a  higher  order  than 
the  differential  operator  occurring  in  the  differential 
equation  that  is  to  be  solved.  This  procedure  introduces 
spurious  particular  solutions.  In  order  for  such  a  scheme 
to  be  stable,  these  spurious  solutions  must  be  damped.  If 
one  uses  third  degree  shape  functions  in  a  manner  analogous 
to  case  2a  to  solve  the  differential  equation  for  the  time 
dependence,  then  one  uses  information  for  u^,  u^,  and 

uk_^.  To  define  an  initial  value  problem  for  the 
differential  equations  the  values  of  uk  and  uk  are  sufficient. 
One  recognizes  that  this  form  of  discretization  will  introduce 
spurious  particular  solutions  and  because  of  the  complete 
symmetry  which  exists  between  the  positive  and  the  negative 
time  direction,  there  exists  for  each  solution  which  is 
damped  another  one  which  is  excited  (that  is,  damped  in  the 
negative  time  direction) . 

Because  of  this  observation,  higher  order  elements  of 
the  character  discussed  in  cases  2  and  3  must  be  rejected. 

For  first  order  elements  this  phenomenon  does  not 
occur,  but  even  then  one  might  have  reservations  because 
of  the  choice  of  the  weight  functions.  In  a  marching 
procedure  one  considers  the  solution  as  known  up  to  a  point 
t^,  and  continues  it  through  the  interval  between  tk  and 
tk+^.  This  step  has  no  influence  on  the  solution  in  the 
preceding  interval.  However,  the  weights  used  in  cases  2 
and  3  take  into  account  the  residuals  between  the  points 
tk_i  and  tk  as  well  as  those  between  the  points  t-k  and 
tk+1«  The  solution  in  the  interval  between  tk  and  tk+1 
is  chosen  in  such  a  manner  that  its  residual  counteracts 
the  effect  of  the  residual  in  the  preceding  interval,  over 
which  one  has  no  control.  One  observes,  on  the  other  hand, 


that  according  to  the  above  stability  analysis  this  has  no 
obvious  undesirable  effect.  For  Courant  number  3  the  methods 
even  give  a  cancellation  of  the  discretization  errors  in 
time  and  space  directions. 

Let  us  now  consider  procedures  for  which  this  objection 
does  not  hold. 

First  we  discuss  shape  functions  which  are 
piecewise  linear  in  space  as  well  as  in  time.  In 
the  choice  of  the  weights  for  the  space  direction  one 
will  observe  the  symmetry  which  exists  between  positive  and 
negative  values  of  y.  Regarding  the  time  direction  a 
corresponding  symmetry  is  not  required.  For  higher  order 
elements  it  is  even  undesirable  (because  of  the  occurrence 
of  spurious  solutions,  as  we  have  seen  above).  Specifically 
we  choose,  for  piecewise  linear  shape  functions 

-£  4  £>°> 


The  weight  function  straddles  the  point  t^  but  not  the  point 
tk+l*  Th;i-S  is  necessary  because  one  will  have  a  jump  of  u 
at  the  point  t.  which  generates  a  5  -  function  when  one  forms 


d^y/dt^.  One  has 
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One  thus  obtains 


a  -  2  u.  ■¥■  u>  +  x  ^  -  o 
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Again,  we  introduce  an  amplification  factor 


(69) 


This  leads  to 
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For  small  values  of  vt^  the  roots  are  conjugate  complex. 
The  absolute  value  can  be  computed  from  the  last  formula. 

It  can  be  obtained  more  simply  by  a  comparison  of  the  first 
and  the  last  term  in  Eg.  (70) .  One  obtains 

-//i 

/f/~  f*  /  (72) 
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A  double  root  is  obtained  for  (vhfc)  =  16.  This  gives 


/ 
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For  (vh^)  >  16  one  obtains  two  real  values  of  p.  The 

particular  solutions  are  always  damped  if  |p|  <  1. 


We  examine  for  which  values  of  vhfc  one  will  have 
|  p |  ®  1.  One  obtains  from  Eg,  (70) 


and 


(vh^)  *  0  ■foi  fm  / 


OO  fry  p*-t 

The  particular  solutions  are  damped  for  all  finite  values  of 
vhfc.  As  vhfc  tends  to  infinity  one  obtains  particular  solutions 
which  alternate  between  positive  and  negative  with  little 
change  of  amplitude.  Fig.  24  shows  exp(iyhfc)  and  a  curve 
p(vhfc)  in  the  complex  plane  for  0  <  yhfc  <  it.  Corresponding 
points  are  connected  with  each  other.  The  useful  range  for 
values  of  yht  is  rather  limited  to  about  0  <  yhfc<  tt/2.  Beyond 
this  point  one  obtains  strong  damping  and  a  considerable 
phase  error.  The  method  can  be  used  even  if  yht  is  large, 
because  such  waves  do  not  create  an  instability,  but  the 
contributions  of  waves  for  large  values  of  yhfc  are  meaningless. 

Now  we  study  again  Eq.  (54)  but  represent  the  unknown 

in  the  interval  from  k  to  k  +  1  by  third  degree  polynomials 

They  are  determined  by  the  known  quantities  u^  and  u^  and 

by  the  unknown  quantities  u.  ,,  and  u’  .  By  this  characteriza- 

^  i  Jc+1 

tion  continuity  of  the  first  derivates  at  the  grid  points  is 
guaranteed.  The  functions  u  represent  the  coefficients  of 
the  eigenfunctions  which  arise  from  a  representation  of  <J> 

(in  the  original  partial  differential  equation)  along  lines 
t  =  constant  in  terms  of  piecewise  third  degree  polynomials 
in  y.  They  are  therefore  expressible  in  terms  of  <f>  and 
at  the  grid  points.  The  continuity  of  the  derivatives  u* 
therefore  implies  continuity  of  <f>t  and  not  only  at  the 
grid  points  but  along  lines  t  =  const.  If 'the  function  <|>t 
is  given  along  a  line  t  =  constant  (which  is  the  case  in  a 
properly  set  initial  value  problem)  then  one  can  determine 
from  it  also  it  does  not  contradict  the  notion  of  a 

properly  formulated  initial  value  problem,  if  we  initially 
assign  the  value  of  u.  If  one  takes  triangular  elements 
and  third-degree  polynomials  (with  the  powers  of  y  and  t 


45 


»■».  it*-  wmfljW 


t 


counted  together)  then  one  has  agreement  of  the  gradient 
in  adjacent  intervals  only  at  the  grid  points,  In  the 
present  formulation  we  have  piecewise  bicubic  elements  and 
then  continuity  of  the  gradient  is  guaranteed  everywhere 
along  the  element  boundaries. 

The  weight  functions  are  to  be  applied  only  to  the 
interval  from  k  to  k  +  1.  No  overlap  over  the  grid  point 
is  needed,  because  of  the  continuity  of  the  first  derivatives. 
In  each  interval  we  have  two  unknowns,  uk+^  and  uk+^.  One 
therefore  needs  two  weight  functions .  We  choose 

w  =  1  for  t^  <  t  <  tk  +  c  hfc  (0  <  c  <  1)  and 


w  =  1  for  tk  +  o  ht  <  t  <  tk+1 


or  equivalently, 

w  =  1  for  tk  <  t  <  tk+1 

w  =  1  for  tk  <  t  <  tk  +  ch 

In  representing  the  space  dependence  of  4>  we  have  also 
considered  examples  in  which  the  weight  functions  were  formed 
by  linear  combinations  of  the  elemental  shape  functions. 
Something  similar  can  be  done  here.  One  could  for  instance 
choose 


W  =  +  N3 


w  =  N2  +  N4 


where  the  functions  are  given  by  Eq.  (6)  with  Ay 
replaced  by  At.  The  first  of  these  weight  functions  is  1. 
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This  possibility  has  not  been  discussed.  The  results 
would  probably  be  of  the  same  character  as  those  for  the 
weight  functions  defined  above.  Let 


At  =  t 


Then,  using  Eq.  (6}  one  has  the  following  representation  for 

u  in  the  interval  between  t,  and  t,  , . 

k  k+1 


-  %  \i-  (■ 


* 


it/ 


•t  /"  <*  *6  (j^lj  *  <\fl  />t  (~i  */(#>] 


Substituting  these  expressions  into  Eq.  (54)  ,  multiplying 
by  the  weight  function  1  and  integrating  over  ^  from  0  to 

.  2  t 

c,  one  obtains  after  multiplication  by  ht 


/  % 

*  ic'k  ht  +  + 

+  *V/  %fc~  c*+  %,  Lt3~  T°y 

*4**  Hr c'-f  + 


ST  O 


(74) 


y 

r;" 
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U,  k 


(74) 

Continued 


+  '1  ' u 


o 


The  second  equation  is  obtained  by  setting  c  =  1.  Let 


Then  one  obtains  the  following  system  of  equations  which 
connects  the  vectors  (uk+1,  uk+1)+  and  (uk»  uk)+. 


L 


'-•*£) 
^  ■/  ^ 


*44 


2. 


c,+H 


iu %  M 
Ml  nt 


=  0 


(75) 


with 


a^  =  -6c  +  6c 

2 

a  2  =  -4c  +  3c 

b1  =  c  -  c3  +  ^c4 

K  _  12  23  .  14 

b2  ■  r3  "3°  +  7° 


c^  =  6c  -  6c 

2 

C2  =  -2c  +  3c 
d!  =  c3  -  jc4 

*2  -  ♦  i*4 


(76) 


The  amplification  ratio  p  defined  by 


48 


•&*t 


is  determined  by  the  following  equation 


0-  ir> 


This  quadratic  equation  for  p  must  be  solved  for  different 
values  of  X  and  c.  The  amplification  ratio  has  either  two 
real  or  two  conjugate  complex  roots.  Eq.  (77)  evaluated  in 
detail  gives 

eft-t){l-6  icO-c))- 

+  [12.  +  A  (~s  -  c(t-c)  +  £-0+  ¥c//-C))J  f  (78) 

,/-6  /  V-f  '  tcO-c»- ^!r]}  '  0 


This  equation  must  be  solved  for  different  values  of  c  and  X. 

The  factor  c(l-c)  is  of  course  unimportant.  One  notices  that  the 
.  2 

coefficient  of  p  and  the  coefficient  of  the  constant  terms 
are  interchanged  if  one  replaces  c  by  (1-c)  and  that  the 
coefficient  of  p  remains  the  same.  Accordingly,  if  c  is 
replaced  by  (1-c),  then  a  solution  of  p  is  replaced  by  p  1. 

For  c  =  j  the  product  of  the  roots  p  is  1.  In  this  case 
the  solution  is  stable  only  if  the  roots  are  conjugate 

complex.  The  discussion  for  c  =  i  best  goes  back  to  the 

.  .  .  ^  1 
original  determinant  Eq.  (77).  Setting  c  =  and  taking  twice 

the  second  row  minus  the  first  row,  one  generates  the  effect 


of  a  weight  function  which,  is  symmetric  with  respect  to 
At/1^  =  i  .  One  obtains 


(~i  *£*)(*+?) 


Hence 


0-  ?)*’  j_  >(i~  Zgl) 

777e?L  "T 


For  A  small  one  obtains 


and  hence 


P  =  1  ±  i<h. 


(as  expected) .  The  right  hand  changes  signs  for  A  *  48, 

A  =  12 ,  and  A=  9.6. 

p  is  complex  for  0  <  A  <  9.6  and  12  <  A  <  48 
p  is  real  for  9.6  <  A  <  12  and  48  <  A. 

The  values  of  (arg  pJA"1^2  versus  A 1/2  for  0  <  A1/2  *  vhfc 
<  9.6  are  shown  in  Fig.  25,  One  sees  that  in  this  stable 
range  the  ideal  value  1  is  fairly  well  approximated.  The 
method  is  unstable  for  c  *  ^  and  P  real. 


One  has,  from  Eq.  (78)  for  general  c  and  two  con 
jugate  complex  roots  of  p 


///* 


t  -  —c(^O)  +  Jx,y  cx 


(80) 


The  first  two  terms  in  the  numerator  and  in  the  denominator 
are  the  same;  for  0  <  c  <  1  they  are  positive.  The  last 
terms  show  that  | p (  <1  for  c  >  •  We  restrict  our 

attention  to  this  region  1/2  <  c  <  1.  Fig.  26  shows  for 
c  =  3/4  the  curves  in  the  complex  p  plane  which  arise  if  X 
is  varied.  They  consist  of  arcs  in  the  complex  plane  which 
are  symmetric  to  the  real  axis.  According  to  Eq.  (80)  these 
arcs  lie  within  the  unit  circle  if  1/2  <  c  <  1.  They  end 
in  branch  points  which  lie  on  the  real  axis  and  are  connected 
by  pieces  of  the  curve  for  which  p  is  real.  These  portions 
may  pass  through  the  unit  circle.  This  can  happen  only  for 
p  =  +  1.  One  obtains  from  Eq.  (76)  for  p  =  1 

£ 

-  c(l-c)  -  6A  =  0. 

2 


Hence 


X  =  0 


and 


(81) 


The  first  value  represents  the  beginning  of  the  curve.  The 
last  value  is  a  second  transition  through  the  point  1.  Only 
for  c  ■  1  will  this  point  lie  at  X  ■  <*>. 
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One  obtains  for  p  =  -1 

-  £-[i  +  tO-o]  *  -z  +  -o 

Hence 

J  *  /i 

and 

J  i+cO'C) 


These  points  coincide  for  c  =  1;  that  is,  the  p-curve  reaches 
the  unit  circle  at  p  =  -1  but  returns , immediately .  The 
arrows  in  Fig.  26  give  the  direction  of  increasing  values 
of  X.  The  curve  starts  with  X  =  0  at  p  =  1  and  moves  along 
two  branches  through  the  complex  plane  to  a  branch  point 
on  the  negative  real  axis.  It  is  continued  along  the  real 
axis  in  the  positive  and  negative  directions.  The  branch 
extending  in  the  negative  direction  passes  through  the  unit 
circle.  Both  branches  double  back  and  meet  at  a  branch  point 
along  the  negative  real  axis.  From  there  the  curve  is  con¬ 
tinued  through  the  complex  plane  in  two  branches  which  meet 
at  a  branch  point  along  the  positive  real  axis.  The  final 
continuations  follow  the  real  axis  in  positive  and  negative 
directions.  One  of  the  branches  passes  through  the  unit 
circle  (at  p  =  1)  .  For  X  »  these  two  branches  end  at 
points  of  the  positive  real  axis.  For  c  =  1  (Fig.  27)  one 
of  the  endpoints  lies  on  the  unit  circle. 

According  to  the  above  analytical  discussions  the  curves 
for  different  values  of  c  have  the  same  character.  All  have 
unstable  regions  (|p|>  1),  except  if  c  =  1.  Stability  for 
c  1  can  be  guaranteed  by  limiting  the  Courant  number 
(X  <  12/(1  +  c(l-c).  For  c  =  1/2  this  limit  is  X  <  12. 


52 


2  2 

Notice  that  these  values  are  close  to  X  =  (vhfc)  =  tt  , 

If  one  imposes  such  a  limitation  then  the  choice  c  -  1/2 
seems  to  be  preferable;  according  to  Fig.  25  it  gives  a  good 
approximation  to  the  propagation  velocity  of  different  waves 

In  comparison  to  linear  shape  functions  the  approxi¬ 
mation  by  third  degree  shape  functions  is  much  better  even 
for  other  values  of  c.  This  is  seen  by  comparison  of 
Figs.  26  and  27  with  Figure  24,  but  as  we  mentioned,  one 
must  limit  the  Courant  number. 

It  is  clear  from  the  outset  that  for  Courant  numbers 
that  are  considerably  above  1  the  solution  of  Eq.  (54)  will 
be  inaccurate.  The  Green's  function,  which  describes  the 
propagation  of  errors,  is  an  oscillating  function.  We  have 
imposed  an  averaging  procedure  for  the  residuals  without 
taking  this  property  of  the  Green's  function  into  account. 
One  is  therefore  resigned  to  the  fact  that  from  a  certain 
value  of  vhfc  on  the  results  will  be  meaningless.  One  hopes, 
however,  that  such  solutions  will  vanish  automatically, 
because  they  are  damped.  The  above  analysis  shows  that  this 
is  not  the  case  for  third  degree  shape  functions  except 
for  c  =  1.  In  this  limiting  case  one  has  as  one  of  the 
weight  functions  constant  weight  throughout  the  interval 
and  as  second  weight  function  a  6  function  at  the  point 
tk+l;  th*s  amounts  to  a  combination  of  weighted  residual 
method  and  a  collocation  method.  Of  course,  a  collocation 
method  has  much  in  common  with  a  finite  difference  approach. 

We  notice  in  passing  that  an  implicit  method  does  not 
automatically  lead  to  a  stable  procedure. 
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SECTION  VII 

A  MORE  COMPLICATED  SAMPLE  PROBLEM 


In  the  preceding  sections  we  have  seen  that  in  many 
methods  stability  can  be  achieved  by  limiting  the  Courant 
number.  There  are  methods  which  are  stable  even  for  very 
large  Courant  numbers,  but  this  entails  a  considerable  loss 
of  accuracy  even  for  fairly  long  waves. 


In  a  problem  with  variable  coefficients  it  is  possible 
to  satisfy  a  Courant  number  limitation  throughout  the  whole 
field  provided  that  the  discriminant  which  determines  whether 
the  problem  is  elliptic  or  hyperbolic  is  bounded  away  from 
zero.  If  the  step  in  the  space  direction  is  fixed,  one 
simply  makes  the  time  step  sufficiently  small.  Of  course, 
for  certain  parts  of  the  field  the  Courant  nuinber  will  then 
be  unnecessarily  small  and  the  progress  in  the  t  direction 
unnecessarily  slow.  At  first  glance  it  seems  that  this  is 
no  longer  possible  if  this  discriminant  vanishes  locally  as 
it  happens  in  the  transonic  problem.  (In  the  present 
discussion  which  is  restricted  to  hyperbolic  equations,  this 
can  happen  only  at  the  boundary  of  the  field.)  We  study  a 
simply  example  of  this  kind. 


Consider  the  partial  differential  equation 

2>v  „ 

-y  *  Jpr’O 

with  boundary  conditions 


(83) 


«  o  fo  r  /  -  O 

<f>  .  O  ^  y  -  / 


(84) 


The  discriminant  vanishes  at  the  boundary  y  =  0  of  the  field. 
In  a  transonic  problem,  x  is  the  downstream  direction.  It 
corresponds  to  t  in  the  previous  examples. 
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For  a  comparison  we  shall  treat  the  wave  equation  also 

_  2V  +  -  o  (85) 

d  x*~ 

with  the  same  boundary  conditions. 


For  simplicity  a  semidiscretized  difference  method  with 
L  +  1  equal  intervals  in  the  y  direction  (L  intermediate 
points)  is  chosen.  Then  one  has 


(86) 


Let 

h<*>  -  H*,  > 


(87) 


The  semidiscretized  finite  difference  approximations  to 
Eqs.  (83)  and  (85)  are,  respectively 


* 

and 


(transonic  equation) 


■  *%  ft*** 


(wave  equation) 


Setting 


one  obtains  the  following  eigenvalue  problems 


with 


and 


j  - 


(transonic  equation) 


*%  +  ~2%  +  -  0 


with  J  <r  >*£** 


(wave  equation) 


(88) 

(89) 

(90) 

(91) 

(92) 

(93) 

(94) 
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The  n  eigenvector  for  the  wave  equation  has  the  elements 


Os 


(*v) 


Sk 


L *  / 


(95) 


Then  one  has 


c*-/ 


A  =  +  Jc'iv ( 


H  F 

L  + 1 


(96) 


The  largest  eigenvalue  is  obtained  for  n  =  L.  If  L  is 
increased,  then  one  obtains  more  eigenvalues  but  within  the 
same  interval  0  <  X^  <4. 


For  the  transonic  problem  one  has  a  generalized  matrix 
eigenvalue  problem 


-2 


l 


l 

-2 

l 


l 


-2 


1-2  1 
1  -2 


1 


- 

—  — ■ 

l 

al 

a2 

3 

a3 

L-l 

1 

a 

1 

L-l 

L  - 

J 

ar. 

(97) 


This  problem  has  been  solved  for  L  =  8,  L  =  12  and  L  =  20 
by  means  of  a  routine  available  in  the  EISPACK  library.  The 
eigenvalues  are  listed  in  Table  1. 

For  the  wave  equation  as  well  as  for  the  transonic 
equation  there  exists  a  maximum  eigenvalue  (X  =  2.4  and  X  =  4, 
respectively) .  For  the  integration  with  respect  to  time  (in 
which  the  question  of  stability  is  decided)  one  must  consider 
the  value  of  v  which  is  connected  to  X  by  Eqs.  (92)  and  (94; , 
respectively.  In  the  preceding  section  we  found  as  a 
stability  limit  for  approximation  of  the  time  dependence  by 
third  order  polynomials 
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TABLE  1 


EIGENVALUES  FOR  DIFFERENT  VALUES  OF  L 
IN  THE  TRANSONIC  CASE 


L  =  8 

L  =  12 

L  =  20 

2.3880 

1.0900 

.70560 

.51691 

.37290 

.22900 

.00672 

.025700 

2.3880 

1.0900 

0.70566 

.52166 

.41367 

.34015 

.27464 

.20644 

.13978 

.081353 

.036374 

.0085800 

2.3880 

1.0900 

.70566 

.52166 

.41376 

.34284 

.29267 

.25531 

.22632 

.20215 

.17847 

.15318 

.12699 

.10112 

.076658 

.054498 

.035403 

.019995 

.0087602 

.0020425 
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TABLE  2 


THE  VALUES  OF  a.  ,  k  =  1...8  (FOR  THE 
LARGEST  EIGENVALUE  OF* THE  TRANSONIC  EQUATION 


D 

L  =  8 

L  =  12 

L  =  20 

1.00000 

1.0000000 

1.000000 

-.38801545 

-.38801545 

-.038801545 

.077142863 

.77142863 

.077142863 

-.1035388 

-.010353876 

-.010353876 

5 

.0010502475 

.0010502475 

.0010502475 

6 

-.000085665 

-.000085665 

-.000085665 

7 

.0000058444 

.0000058449 

.0000058449 

8 

-.000000341697 

-.00000034273133 

-.00000034273133 

TABLE  3 


THE  VALUES  OF  ak ,  k  =  1 ...  8 ,  FOR  THE  SECOND 
LARGEST  EIGENVALUE  OF  THE  TRANSONIC  EQUATION 


w 

CO 

II 

L  =  12 

L  =  20 

1 

-.85928478 

-.85928475 

-.85928475 

2 

-.78197087 

-.78197083 

-.78197083 

3 

1.000000 

1.0000000 

1.000000 

4 

-.48795330 

-.48795340 

-.48795340 

5 

.15152046 

.1512073 

.15152073 

6 

-.034773134 

-.034773974 

-.034773974 

7 

.0063442942 

.0063478456 

-.0063478456 

8 

-.0009411981 

-.00096327439 

-.00096327439 

TABLE  4 

THE  VALUES  OF  a.  ,  k  =  1 ...  8  FOR  THE  FOURTH 
LARGEST  EIGENVALUE 


w 

L  =  8 

L  =  12 

L  =  20 

1 

.67082316 

.67489092 

.67489121 

2 

.99489036 

.99772005 

.99772028 

3 

.29041777 

.27961324 

.27961259 

4 

-.86441535 

-.87608045 

-.87608122 

5 

-.2319448 

-.20371902 

-.20371726 

6 

1.0000 

1.0000000 

1.0000 

7 

-.86952195 

-.9262249 

-.92622843 

8 

.40721511 

_ 

.52975421 

_ 

.52976194 

To  improve  the  accuracy  one  may  choose  a  somewhat  lower  limit. 
Substituting  Eqs.  (92)  and  (94)  into  this  condition  and  using 
the  respective  maximum  values  of  X  one  obtains 

4*-  (98) 

(transonic  equation) 

(99) 

4**  ^.*4*  ,  . 

*  J  (wave  equation) 

Important  are  the  powers  of  h^  (3  for  the  transonic  equation, 

2  for  the  wave  equation) .  One  must,  of  course,  reduce  h 
as  one  reduces  h^,  but  more  strongly  for  the  transonic  4 
problem  than  for  the  wave  equation.  But  even  here,  a  mesh 
for  which  the  time  integration  converges  can  always  be  found. 

One  observes  that  in  the  transonic  example  the  largest 
eigenvalues  for  different  values  of  L  agree  nearly  perfectly. 
This  does  not  happen  for  the  wave  equation.  Also  the 

/„  v 

components  a^  of  the  eigenvectors  (normalized  by  the  require¬ 
ment  that  the  largest  component  of  the  vector  with  components 
(n) 

a^  '  be  one)  agree  surprisingly  well  for  different  values 
of  L.  (See  Table  2  for  the  largest  eigenvalue, 

Table  3  for  the  second  largest  and  Table  4  for  the  fourth 
largest.)  For  large  values  of  k  the  components  a^  of  these 
vectors  become  extremely  small,  therefore  only  the  components 
for  k  =  1  to  k  =  8  are  included  in  these  tables.  One  sees 
that  the  effect  of  the  eigenvectors  pertaining  to  the  largest 
eigenfunctions  is  restricted  to  the  immediate  vicinity  of  the 
parabolic  line.  As  L  is  increased  and  the  interval  h^  becomes 
smaller,  the  region  where  the  largest  eigenvalues  are  of 
importance  contracts  because  close  agreement  of  the  eigenvectors 
occurs  if  one  compares  the  same  values  of  k  (not  of  y  =  k  hy) 

This  analysis  shows  that  the  Courant  number  limitations 
are  not  as  serious  as  one  might  assume  at  a  first  glance. 

If  the  contributions  belonging  to  the  largest  eigenvalues 


should  become  unstable,  then  this  instability  will  be 
confined  to  the  immediate  vicinity  of  the  parabolic  line 
y  =  0. 

The  eigenvectors  have  been  found  by  direct  computation. 
The  is  probably  the  most  practical  way  even  if  one  does  not 
use  the  EISPACK  routine,  which  is  somewhat  oversophisticated 
for  the  present  simple  problem.  The  eigenfunctions  can  be 
expressed  in  terms  of  Bessel  functions.  This  may  be  useful 
if  one  wants  to  understand  their  asymptotic  properties,  for 
instance,  the  fact  that  the  largest  eigenvalues  are  nearly 
identical. 

The  system  Eq.  (97)  can  be  regarded  as  a  three  point 
recurrence  relation 

4  4  o 

It  resembles  the  recurrence  relation  for  Bessel  functions 

where  Z  represents  a  linear  combination  of  J  and  N  with 
P  P  P 

coefficients  that  are  independent  of  p.  We  set 

Then  one  has 

o  (100) 

To  identify  the  two  relations  one  must  postulate 

4^  -2  -  i* 

x 

The  order  p  must  increase  by  1  as  k  is  increased  by  1. 
Therefore, 

^ 

and,  using  this  result 

/6  ■  4  ~  2/J 

Thus,  one  finds  that 

4  - 


(101) 


satisfies  the  recurrence  relation  Eq.  (100) .  The  eigenvalue 
X  is  determined  by  the  boundary  conditions 

Accordingly,  one  must  determine  X  and  the  linear  combination 

between  J  and  N  in  such  a  manner  that 
P  P 

£  0  (102) 

and 

£  Ufa)*  O  (103) 

If  L  is  very  large,  then  the  order  L  =  1  — ( 2/X )  is  very  large 
in  comparison  to  the  argument  (2/X)  of  the  Bessel  functions. 
Bessel's  equation  reads 

If  the  order  is  very  large  in  comparison  to  the  argument, 
then  (except  for  a  constant)  Zp  =  X-p.  It  then  follows 
from  Eq.  (103)  that  Z  is,  in  essence,  given  by  Jp  =  const.  Xp 
while  the  contribution  of  N  is  very  small  (except  for  k  -  L) . 
This  explains  why  the  eigenfunctions  are  so  close  to  each 
other  for  small  values  of  k. 
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SECTION  VIII 


A  RATIONALE  FOR  THE  CHOICE  OF  THE  WEIGHT  FUNCTIONS 


So  far  we  have  chosen  the  weight  functions  on  intuitive 
grounds  and  then  examined  their  effect  on  the  .stability. 

In  this  section  the  attempt  is  made  to  provide  for  a  more 


rational  basis. 


We  ask  for  the  best  possible  choice  of  a  weight  function 
for  the  differential  equation  Eq.  (54)  assuming  that  v  has 
a  fixed  value.  As  usual,  we  replace  u  by  some  approximation 
using  shape  functions  which  are  related  to  certain  data 
within  the  flow  field;  in  our  case  these  are  the  values  of 
u  and  u‘  at  the  grid  points.  The  approximation  so  obtained 
is  denoted  by  u.  With  weight  functions  w.^  and  w2 ,  one 
obtains  the  following  equations ,  from  which  the  values  of 


u^  and  u£  are  determined. 


/  fa  t)  -  o, 


(104) 


The  specific  choice,  to  be  justified  presently,  is 


V,  =  CM(lAt)  .  %  -  f&H'fat  ) 


(105) 


Then  one  has  the  condition 


/'w-  *  * 


Aih,  fa 


Jfc  m  O 


(106) 


Let  Au  be  the  difference  between  the  exact  solution  u  and 


the  approximation  u 


U,  m  it,  +  AU 


One  obtains,  by  substituting  this  expression  into  Eq.  (54) 


with 


1 
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Hence,  provided  that  for  Au  =  0,  Au*  =  0  for  t  =  t^ 

t  * 

v'[~Cd(VAt) J" &(t)'U*v(V4T) tfc  +  J  fc  ft )&$(»&&&] 


K 


6  iz 

4  1 


(107) 


Hence,  with  Eqs.  (106) 

s  °>  J-  ° 

One  sees  that  with  the  choice  Eq.  (105)  of  the  weight  function 
the  errors  in  u  and  u‘  vanish  at  the  grid  point  t^+^,  and  by 
induction  at  all  subsequent  grids  points. 

The  argument  is  slightly  different  for  the  space 
dependence.  The  starting  point  is  the  equation 

#L 


J>£  2# 

Ux  lyK  * 


(108) 


We  compare  an  exact  solution  and  an  approximation  which  have 
the  same  dependence  upon  time.  The  exact  solution  is  given 
by 

?  fy,t)  *  Ufy)  tx/i (+  cvt )  «  4*f>(Lyy) tty  fiivt) 

fyCy,t)  me/Ur/JycxtfJcvt;  *  C\>  *Kf>  fry)  **/>(*  ivt) 

One  has  specifically 

-  £*/>(cvkh  ) 

*  6  ’ey 

The  factor  h  1  in  u^1  occurs  because  by  our  original  definition, 
the  prime  denotes  the  derivative  with  respect  to  t/h.  The 
function  u  satisfies  the  equation 

Jh/Jt1-  +  v'to  »  O 

The  approximation  is  assumed  to  have  the  form 

tl(y)  tX/>  Civ  t) 


64 


Now  it  is  assumed  that  u(y)  is  given  by  piecewise  third 
degree  polynomials  with  continuity  of  the  function  and  of 
the  first  derivative.  Let 


and 


-  i(kh) 


With  tentative  weight  functions  w^(Ay)  one  finds  as  conditions 
from  which  u^  and  u^'  are  determined 

/  '  0  (109> 

Previously  we  proceeded  as  follows.  Particular  solutions 
are  obtained  by  setting 


a  C  ) 

U-fJ  *  CC  (iMskh ) 

This  leads  to  a  homogeneous  system  for  c  and  C.  The  elements 
of  the  governing  matrix  depends  upon  y.  The  vanishing  of  its 
determinant  establishes  a  relation  between  y  and  v  (no  matter 
how  the  weight  functions  are  chosen) .  We  set 

A  H'fjl)  *  -  £(y) 

One  then  obtains 


A  ic  +  u,  +  & 


f 

ey*- 

A  0  Co  %  +* 

4  ~  ~Zji  *  p  ^ 

Then  y  , 

Jmh 

My)jt*,(»Ay)</y  -  &»fay)f/l(y) e*  fa yj/y/ 
J5t  j 

alc  -Ait  *  - Jty»(»cy  J fc(u)jih,(»ay)yy  - cas{ vAy) f  k(y) cat  (My) ^y 
4  /A 


(110) 
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Now  we  impose  the  condition  that  the  errors  in  Au  and  Au' 
vanish  at  the  grid  points  in  the  y-direction  the  left  sides 
of  the  last  equation  vanish.  In  order  for  the  right  side  to 
vanish,  it  is  necessary  that  the  integrals  vanish  separately. 


+  f&)ces()>4  y)Jy 


d 


These  conditions  are  identical  with  Eqs.  (109)  if  one  chooses 
«=• 

(111) 

4^  -  dosCvA!/) 

If  these  conditions  are  satisfied,  one  obtains  because  of 
Eq.  (108) 


and  then  it  follows  that  u(y) 
It  does  not  matter  how  u 
provided,  of  course,  that  if 


is  periodic  iwy  with  period 
(y)  is  related  to  uk  and  u£ 


then  „  . 

%(***) 

The  initial  data  will,  in  reality,  contain  linear 
combinations  of  waves  exp (ivy)  with  different  values  of  v. 

The  weight  functions  sin(vAy)  and  cos  (vAy)  are  perfect  only 
for  one  of  them?  for  all  other  wave  lengths  some  error  will 
occur.  At  befit,  it  is  possible  to  tune  the  weight  functions  to 
a  certain  wave  length,  then  it  will  be  nearly  correct  for 
adjacent  frequencies.  The  weight  function  for  very  long 
waves  are  given  by 
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and 


li**,  Cfii  (v  a  y  )  •  / 

V-*  c> 

/t+v  (v'jvi/fvag)  -  ay 
\>-*6 

Figure  28  is  similar  to  Figures  12  and  13.  It  shows 
for  third  degree  weight  functions  the  wave  speed  for 
different  wave  lengths  with  this  choice  of  the  weight  functions. 
For  short  waves  (large  values  of  y  or  v)  the  result  is 
certainly  not  better  than  the  weight  functions  on  which 
Figures  12  and  13  are  based.  It  might  be  preferable  to 
attune  the  weight  function  to  a  value  of  v  somewhat  larger 
than  zero . 

The  Duhamel  solution  of  the  inhomogeneous  problem 
used  in  Eqs.  (107)  can  be  regarded  as  a  representation  of 
the  Green's  function  for  the  ordinary  differential  equation  (54) 

The  Green's  function  for  the  wave  equation  can  be 
obtained  by  a  Fourier  analysis  with  respect  to  the  y  and  t 
directions.  The  circular  frequency  v  will  then  vary  from 
0  to  infinity.  The  weight  functions  which  one  would  obtain 
by  applying  the  above  rationale  are  w^  =  1,  Wj  =  At,  w3  =  Ay 
and  w.  =  Ay  At.  They  give  equations  for  <f>,  <j»  .  <j>.  and  <)>  . 

4i  y  t  yt 

at  each  grid  point.  These  weight  functions  are  attuned  to 
the  low  freqCiency  components  of  the  Fourier  decompositon. 

The  Green's  function  for  the  Laplace  equation  is 
given  (except  for  a  constant)  by  log(r-r').  One  then  obtains 
as  the  effect  of  a  residual  R(r') 

iff)  -  Jf  *(?'!  iylr-T'l  A‘Jf‘ 

From  this  expression  one  can  derive  weight  functions  by  the 
requirement  that  the  long  distance  effect  of  the  residual 
from  each  element  be  small.  Let  the  element  center  be  at 
r'  =  r  with  coordinates  x  =  xQ,  y  =  yQ.  Then, 
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|r'  -  rQ|is  small  in  comparison  toCr  -  rQ|.  One  obtains,  by  a 
development  of  log(r'  -  "r)  with  respect  to  (x'  -  x  )  and 

<y  - 

ff  A(r')  Jy  (r -r' /</*'(*','  -  (?')c/x'c/y' 

//urn*'- VMS 

/  /tffker  a  re/e  r  term: 


This  suggests  three  weight  functions  for  each  element,  namely, 

W1  =  1,  w2  “  A*  and  W3  =  Ay.  (One  remembers  that  for  the 

Laplace  equation  and  bicubic  elements  one  expresses  the 

solution  in  terms  of  <f>,  <f>  and  <J>  for  each  grid  point.) 

x  y 

The  number  of  weight  functions  obtained  by  the  above  reasoning 
is  correct.  There  is,  of  course,  a  question  whether  one  should 
use  these  weight  functions  instead  of  those  which  arise  from 
an  extremum  formulation.  The  same  approach  can  be  used  to 
obtain  weight  functions  for  the  three  dimensional  Laplace 
equation. 


In  a  corresponding  manner  one  can  discuss  the  choice 
of  the  weight  functions  for  the  three  dimensional  Helmholtz 
equation.  There  one  has  as  effect  of  a  residual  (except  for 
a  constant) 


Jf/  X'(x'ii/>jZ’) 


/r*-  *7 


The  weight  functions  are  obtained  from  an  approximate 
representation  of  the  Green's  function 


CjriM'l?-  r'/  Cnfa/r-r0f)  Cnju>/r-rt/) r  , 

Tf^F  /«-*#' -V  **-*&-»* 

-  ~  — [(*-**)(*' -*o)  t  fy 

f  fttyker  enrder  terms 
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The  expression  appearing  in  the  second  line  dies  out  only 
as  l/(r-r')»  while  that  in  the  first  line  dies  out  as  l/(r- 
(Either  x  -  xQ  or  y  -  yQ  is  of  the  same  order  as  r.)  The 
above  form  suggests  weight  functions  w^^  =  1,  w2  =  Ax, 
w^  =  Ay,  w^  =  Az,  but  strictly  speaking  the  goal  of  suppressing 
long  range  effects  is  achieved  only  if  u(x'-Xq),  Vi(y'-y0)  and 
u(z'-Zq)  are  small  within  the  element  under  consideration. 

The  size  of  admissible  finite  volumes  must  be  smaller  than 
the  wave  length. 

The  idea  of  making  long  range  effects  small  fails  for 
hyperbolic  equations.  The  effect  of  a  residual  is  given  by 

JJf  WtWM-tr- (r-s'f- 

The  square  root  vanishes  at  the  characteristic  cone  through 
the  point  t',  y',  z',  and  for  this  vicinity  a  development 
of  this  expression  is  no  longer  feasible.  It  is  impossible 
to  supress  long  range  effects.  Some  further  insight  is 
obtained  if  one  considers  the  finite  element  concept  in 
conjunction  with  the  idea  of  characteristics. 
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SECTION  IX 

FINITE  ELEMENT  AND  CHARACTERISTIC  COORDINATES 

The  presence  of  characteristics  imposes  a  d&finite 
structure  to  hyperbolic  equations;  one  expects  that 
this  fact  has  some  bearing  on  the  implementation  of  the 
idea  of  weighted  residuals.  Perturbations  (for  instance  the 
effects  of  a  residual)  propagate  along  characteristics. 

In  Appendix  VI  the  Green's  function  belonging  to  a  two 
dimensional  problem  has  been  studied.  A  local  source  in  a 
flow  with  a  Mach  number  /T  gives  a  perturbation  in  the  velocity 
only  along  the  characteristics  emanating  from  the  point 
at  which  the  source  is  located.  Along  these  characteristics 
one  has  a  step  in  the  potential.  The  perturbation  velocity 
has  the  direction  normal  to  them;  it  is  given  by  a  6  function. 
Also  the  perturbation  in  the  mass  flow  vector  is  given  by 
a  delta  function;  it  has  the  direction  of  the  characteristics. 
The  effect  of  such  a  source  does  not  die  out  with  distance. 

In  choosing  weight  functions  for  elliptic  equations,  one  can 
use  a  plausibility  argument.  The  long  range  effect  of 
truncation  errors  can  be  reduced  by  postulating  that  within 
small  regions  the  residuals  counteract  each  other.  This 
argument  cannot  be  applied  here:  the  effect  of  two  sources 
which  lie  on  different  characteristics  will  not  die  out  with 
distance  (although  it  is  confined  to  a  narrow  region) .  We 
study  here  to  which  extent  this  state  of  affairs  can  be 
taken  into  account  by  the  choice  of  the  element  shape  and 
of  the  weight  functions. 

Consider  plane  and  axisymmetric  flows  at  a  Mach  number 
/2\  Then  one  has,  respectively 


4s"  ° 

(112) 

(113) 
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If  small  values  of  r  are  excluded,  then  it  is  practical  to 
introduce  in  the  axisymmetric  case 


Then  one  obtains 
■>/ 


-  2*  *  * 
0x* 


Introducing  characteristic  coordinates 

t  -  *  -y 


or 


X/r 
r  s  *‘K 

one  obtains 

ay 


and 


in  the  two  dimensional  case. 


or 


4-^ 

Ki-\> 

.  *  -i 


in  the  axisymmetric  case 


(115) 


(116) 


(117) 


(118) 

(319) 


The  axisymmetric  case  has  been  included  because  simplicity 
of  the  plane  problem  might  lead  to  faulty  generalizations. 

Now  we  consider  the  equation 

(120) 

where  f(£,ri)  is  either  considered  as  known  or  given  as  a 
function  of  <)>  and  its  lower  derivatives. 

Assume  that  the  element  boundaries  are  characteristic 
lines.  This  means  that  we  deal  with  a  rectangular  grid  in 
the  £,n  plane.  The  solution  in  one  quadrangel  element  is 
completely  determined  if  one  knows  the  values  of  <J>  along 
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two  adjacent  sides.  From  the  characteristic  conditions  one 
then  obtains  the  values  of  ^  and  <f>  along  the  other  sides 
(<|>  along  the  side  for  which  £  =  const,  cf> ^  along  the  side 
for  which  n  =  const) .  The  potential  is  obtained  by 
integrations.  The  effect  of  the  residual  can  be  judged  by 
considering  the  falsification  of  <J>  (or  of  <{>^and  <j>r)  along  the 
newly  computed  sides  of  the  quadrangle.  This  criterion  is 
preferable  to  an  evaluation  of  the  residual  at  infinity 
because  the  errors  do  not  die  out.  If  the  function  f(£,n) 
is  known,  then  one  has  by  a  direct  integration, 

A 

Jo 

% 

f((L  l’%>  -4  (b  W  *  0 

?* 


The  values  of  ^(C  =  £-^,11)  and  <f>^{£;,ri  =  n1)  (and  the  values 
of  $  which  arise  from  them  by  an  integration)  is  the  only 
information  needed  in  order  to  continue  the  computation 
in  adjacent  elements.  Actually,  one  must  approximate  <f>  by 
means  of  a  finite  number  of  parameters.  Equations  for  these 
parameters  are  obtained  by  applying  weight  functions.  But 
since  only  the  data  along  the  element  boundaries  £  =  and 
n  =  rij  are  needed  it  is  sufficient  to  use  weight  functions 
wi  (n)  along  the  segment  £  =  ^  =  const  and  Wj(£)  along  the 
segment  n  =  =  const.  One  therefore  obtains 

/  *■  ti>K  a-i. iJ  -4  Q-b,  V  -/ AS vWi  -  « 

l6  J  1  A 


and  similarly  for  the  second  equation. 
Hence,  in  an  obvious  manner 


0 


(121) 
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This  is  a  weighted  residual  expression  (The  integration  is 
extended  over  the  whole  element) .  The  weight  Wj  depends  only 
upon  n*  In  the  other  case  one  obtains  weights  which  depend 
only  upon  £• 

In  the  ax i symmetric  cases  we  derived  two  forms  of  the 
differential  equation  depending  upon  which  form  one  chosses 
one  has  as  weight  functions  either 


<  4 ) 

or  V;  {£} 

or  %  ^ 

It  is  assumed  that  the  residual  is  formed  for  the  original 
differential  equation  (113)  in  each  case.  There  is  still 
some  arbitrariness  in  the  choice  of  the  weights.  For  large 
values  of  r  the  difference  is  unessential.  One  has  very 
little  variation  in  one  of  the  characteristic  directions 
while  in  the  other  direction  the  choice  of  the  weight  function 
depends  upon  what  kind  of  details  one  is  willing  to 
disregard. 

In  principle,  the  errors  so  admitted  do  not  die  out 
with  distance.  The  reason  for  disregarding  them  is  that  they 
express  unessential  details.  A  finite  element  representation 
depends,  as  always,  on  a  certain  smoothness  of  the  function 
that  is  to  be  represented.  Discontinuities  of  the  derivatives 
are  admissible  along  the  characteristics.  (They  might  be 
introduced  by  discontinuities  in  the  boundary  conditions.) 

Such  discontinuities  at  the  characteristic  element  boundaries 
are  compatible  with  the  present  formulation  and  if  they 
are  present,  a  finite  element  procedure  based  on  the  elements 
bounded  by  characteristics  may  give  better  results  than  the 
finite  element  procedures  discussed  in  the  preceding  sections. 
To  see  which  steps  one  has  to  take  we  consider  the  differential 
equation  <|>£  +  f(£,n)  =  0  and  assume  that  in  a  characteristic 

quadrangle  the  function  $  is  represented  by  bicubics.  This 
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requires  that  one  knows  at  each  of  the  four  corners  tp ,  <p^, 

<p^  and  <j>£r)  is  given  by  the  differential  equations. 

Known  are  the  data  along  the  lines  1,  2,  and  1,  3  (See  Figure 
29)  . 

Specifically,  <J>1#  <p^ ,  1 ,  >r ,  ± ,  4>2,  *5  2’  ^3'  ^,3 

unknown  are  <|>  2 ,  ^£,3'  ^5  4'  4'  and  ^4 

Accordingly,  one  can  use  five  weight  functions.  One  of  them 
is  obviously  the  constant,  which  depends  neither  on  5  or  n, 
for  the  other  weight  functions  one  will  choose  two  which  depend 
solely  upon  £,  and  two  which  depend  solely  upon  n.  Further 
details  are  omitted.  Some  modifications  are  required  along 
a  line  for  which  initial  conditions  are  given.  This  line  must 
coincide  with  a  characteristic,  and  one  deals  with  triangular 
elements . 

The  stability  analysis  in  all  previous  discussion  has 
been  carried  out  for  the  wave  equation.  In  the  present  case 
one  then  considers  particular  solutions  of 


which  are  given  by 

f  *  6Q)  +  fav 

For  bicubic  shape  functions  one  has  as  representation  in  one 
element 

i>  = 

where  the  g^'s  and  h^ ' s  are  polynomials  of  the  third  degree. 
Among  expression  of  this  form,  there  are  some  of  the  form 
f^(C)  and  f2  (n)  •  This  means  that  the  differential  equation 
will  be  exactly  satisfied.  It  follows  that  the  method  is 
stable  in  the  same  sense  as  one  speaks  of  stability  for  the 
methods  examined  in  preceding  sections.  Here  we  have  neither 
damping  nor  dispersive  errors;  the  only  inaccuracies  arise  in 
satisfying  the  boundary  conditions. 
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Of  course  the  characteristics  are  not  fixed  a  priori 
in  nonlinear  problems.  Accordingly,  one  loses  control  over 
the  shape  of  the  elements.  This  will  probably  make  programming 
more  cumbersome. 

An  analogous  approach  to  the  three-dimensional  problem 
leads  to  serious  difficulties.  In  three  dimensions  one  has 
characteristic  surfaces  which  can  be  chosen  in  a  great  variety 
of  ways.  To  generate  such  surfaces,  one  can,  for  instance, 
choose  in  an  initial  value  plane,  a  two-dimensional  mesh 
consisting  of  quadrangles  and  then  construct  the  characteristic 
surfaces  which  emanate  from  these  mesh  lines.  One  obtains 
two  surfaces  for  each  line.  If  the  lines  are  straight,  one 
would  obtain  two  characteristic  planes  for  each  line.  By 
forming  the  lines  of  intersection  of  these  planes  with  a 
noncharacteristic  plane  nearly  perpendicular  to  them,  one 
obtains  a  pattern  similar  to  those  for  the  two-dimensional 
characteristic  method.  But  one  has  two  families  of  lines 
in  the  initial  value  plane  and  thus,  one  obtains  four 
different  families  of  characteristic  surfaces.  The  finite 
volumes  bounded  by  such  surfaces  are  not  well  suited  for  a 
computation.  Moreover,  these  surfaces  depend  upon  the  choice 
of  the  original  grid  and  a  great  number  of  choices  are 
possible. 

The  characteristics  in  two  dimensions  and  in  three 
dimensions  are  lines  or  surfaces  along  which  discontinuities 
may  propagate.  (Such  discontinuities  are  introduced  either 
by  the  initial  or  by  the  boundary  conditions.)  In  the 
two-dimensional  problem,  characteristics  along  which 
singularities  (or  also  steep  gradients)  propagate  are  readily 
found,  they  are  members  of  the  family  of  characteristics  by 
means  of  which  the  computations  are  carried  out.  In  the 
three-dimensional  case  it  would  be  necessary  to  orient 
characteristic  surfaces  according  to  the  discontinuities  (and 
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steep  gradients)  as  they  are  introduced  by  the  initial  or 
boundary  conditions.  This  may  make  it  necessary  to  change 
the  definition  of  the  elements  as  one  goes  along. 

In  a  typical  three-dimensional  finite  volume  approach 
one  might  consider  it  as  desirable  to  build  up  the  whole 
flow  field  from  finite  volume  elements  which  are  bounded  by 
characteristic  surfaces  in  a  manner  which  takes  the  regions 
of  dependence  into  account.  This  is  very  difficult,  if  not 
impossible.  Assume  for  instance  that  one  has  a  rectangular 
grid  in  a  noncharacteristic  initial  plane.  The  characteristic 
surfaces  described  above  form  four-sided  pyramids  which  lie 
wholly  in  the  region  of  influence  of  the  data  assigned  within 
an  initial  quadrangle.  However,  one  needs  very  complicated 
elements  in  order  to  fill  during  the  next  step  the  spaces 
between  these  pyramids. 

A  finite  difference  formulation  based  on  the  idea  of 
characteristic  surfaces  is  better  able  to  cope  with  this 
problem  because  it  defines  the  approximation  only  at  the 
grid  points. 

In  the  two-dimensional  problem  a  finite  element  procedure 
combined  with  the  method  of  characteristic  can  also  be 
motivated  by  considering  the  Green's  function  (see  Section  ix)  . 
The  Green's  function  of  the  three-dimensional  problem  has  a 
very  complicated  singularity  along  the  characteristic  cone. 

This  makes  it  difficult,  if  not  impossible,  to  derive  from 
it  ,  forms  of  the  weight  functions  which  are  well  fitted  to 
the  problem. 

The  singularity  becomes  more  managable  (in  fact,  it 
becomes  rather  close  to  the  Green's  functions  for  the  two- 
dimensional  problem),  if  along  certain  lines,  the  residuals 
are  smooth  functions  and  one  carries  out  an  integration  over 
these  lines. 
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One  feasible  method  of  approach  can  be  described  as 
follows.  One  introduces  in  an  initial  plane  one  family  of 
nonintersecting  curves.  In  a  plane  problem  these  would  be 
straight  lines  in  the  direction  in  which  <fi  and  the  velocities 
do  not  vary.  In  an  axisymmetric  problem  they  would  be  circles 
in  a  plane  normal  to  the  axis  of  symmetry.  Then  one  forms 
characteristic  surfaces  which  emanate  from  these  lines.  To 
be  specif ic, at  each  point  of  one  of  these  lines  one  constructs 
the  characteristic  cone.  The  characteristic  surfaces  mentioned 
above  are  the  envelopes  of  these  cones.  The  lines  along  which 
the  characteristic  surfaces  are  tangent  to  the  cones  are 
called  bi-characteristics.  One  thus  obtains  elongated 
subvolumes  with  triangular  or  quadrangular  cross  sections 
in  the  initial  elements, and  quadrangular  cross  sections  in 
subsequent  elements.  These  cross  sections  lie  approximately 
in  planes  determined  by  the  bi-characteristics.  In  the  plane 
case,  these  subvolumes  are  rods  with  appropriate  triangular 
or  quadrangular  cross  sections,  in  the  axisymmetric  case  they 
are  rings.  The  characteristic  conditions  contain  only 
derivatives  within  a  characteristic  surface,  one  of  the 
directions  in  which  derivatives  are  formed  is  that  of  the  bi¬ 
characteristic,  the  other  one  in  the  lengthwise  direction  of 
the  rods  or  rings  mentioned  above . 

The  elements  to  be  used  in  a  finite  element  approach 
are  cut  out  from  these  subvolumes  by  planes  approximately 
determined  by  bi-characteristics.  In  the  plane  and  axisymmetric 
case  this  leads  directly  to  the  familiar  characteristic 
conditions.  More  general  problems  can  probably  be  treated 
too,  provided  that  the  characteristic  surfaces  are  oriented 
so  that  they  coincide  with  surfaces  where  singularities 
of  strong  gradients  occur.  If  this  condition  is  not 
satisfied,  then  singularities  or  strong  gradients  will 
probably  express  themselves  in  the  same  manner  as  in  the 


usual  finite  difference  methods,  namely  by  oscillations 
around  the  correct  values.  If  strong  gradients  are  detected 
then  it  is  probably  best  to  make  a  transition  to  a  different 
more  suitably  oriented  set  of  characteristic  surfaces. 

The  author  has  not  been  able  to  study  the  practical 
aspects  of  such  a  procedure. 


SECTION  X 


GENERAL  OBSERVATIONS  AND  SUMMARY  OF 
SPECIFIC  RESULTS 


It  was  already  mentioned  in  the  introduction  that  the 
attempt  to  derive  a  finite  element  method  for  hyperbolic 
problems  from  a  variational  formulation  is  not  likely  to  be 
useful;  at  best  the  idea  is  unnecessarily  confining.  In 
this  report  the  method  of  weighted  residuals  is  advocated 


instead. 


We  add  that  the  idea  of  gaining  an  extremum  formulation 
by  postulating  that  the  sum  of  the  squares  of  the  residual 
be  minimized  is  of  doubtful  value.  Interpreting  this 
method  in  terms  of  weighted  residuals,  one  finds  that  the 
weight  is  proportional  to  the  residual.  But  the  long 
distance  effect  of  residuals  is  linear,  no  matter  how  large 
or  small  they  are.  Therefore,  the  weight  should  be  independent 
of  the  magnitude  of  the  residual.  Minimizing  the  sum  of  the 
squares  of  the  residuals  suppresses  short  waves  (because  they 
give  larger  residuals  at  equal  amplitude) ,  but  the  method  is 
fairly  insensitive  to  long  wave  errors. 


Hyperbolic  problems  differ  from  elliptic  problems  by  the 
fact  that  the*  solutions  do  not  automatically  smooth  out. 

In  an  elliptic  problem  the  waviness  of  the  solution  at  some 
distance  from  the  boundary  reflects  the  local  waviness  of 
the  inhomogeneous  term.  In  a  hyperbolic  problem  the  waviness 
is  determined  by  the  initial  conditions  the  boundary  conditions 
and  the  inhomogeneous  part  in  all  of  the  characteristic 
forecone  pertaining  to  the  point  under  consideration.  Therefore, 
one  cannot  count  on  the  smoothness  of  the  solution. 


One  observes  that  solutions  which  are  wavy  in  the  space 
direction  are  also  wavy  in  the  time  direction.  Even  short 
waves  will  not  be  damped.  This  fact  is  directly  connected 
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with  the  existence  of  characteristic  surfaces.  Along  such 
surfaces  singularities  may  propagate;  their  attenuation 
with  distance  is  only  small.  Assume  that  by  the  initial 
conditions  a  discontinuity  (in  the  second  derivatives,  say) 
is  introduced  in  the  space  direction.  In  order  to  represent 
such  initial  data  by  a  Fourier  analysis  one  needs  all  wave 
length  including  the  very  shortest  ones.  At  a  later  time 
this  singularity  will  still  be  present,  although  in  a 
different  position.  If  short  waves  were  more  strongly 
damped  than  long  ones  such  singularities  would  vanish. 

This  explains  why  stability  is  particularly  critical 
in  hyperbolic  problems.  In  essence,  the  individual 
particular  solutions  have  neutral  stability  (neither  damped 
nor  excited) .  The  error  introduced  by  an  approximation  may 
change  neutral  stability  into  instability.  In  elliptic 
problems  where  the  particular  solutions  die  out  with  distance, 
the  same  kind  of  error  will  falsify  the  rate  of  attenuation, 
but  it  is  unlikely  that  an  exact  solution  which  is  stable 
will  change  into  one  that  is  unstable.  The  situation  is 
particularly  favorable  because  the  contributions  of  those 
particular  solutions  which  are  strongly  falsified  are  small 
to  begin  with  and  furthermore,  because  they  are  most  strongly 
damped . 

The  finite  difference  as  well  as  the  finite  element 
method  presupposes  that  the  solutions  are  smooth  enough  so 
that  they  can  be  approximated  within  the  inidividual  elements 
by  simple  standardized  expressions.  In  elliptic,  as  in 
hyperbolic  problems,  one  must  choose  a  grid  which  is  suitable 
to  represent  the  essence  of  the  desired  solution  but  disregards 
short  wave  roughness.  In  elliptic  equations  the  mesh  size 
satisfying  this  requirement  is  mainly  determined  by  the 
boundary  data.  In  principle,  one  might  use  a  larger  grid 
at  some  distance  from  the  boundaries.  In  hyperbolic  problems 
the  same  waviness  is  to  be  expected  throughout  the  field. 
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The  discussions  carried  out  in  this  report  show  the 
effect  of  truncation  errors.  Let  us  confine  our  attention 
to  the  semi-discretized  methods  investigated  in  Sections  II 
through  VII.  If  the  wave  speeds  a.re  reproduced 
with  sufficient  accuracy  for  the  range  of  wave  lengths 
which  is  important,  then  the  truncation  error  remains  within 
acceptable  bounds.  Shorter  waves  in  the  initial  data  and  in 
the  inhomogeneous  part  are  falsified  because  of  truncations 
errors,  but  it  is  assumed  that  these  components  are  initially 
small.  They  will  remain  small  unless  the  method  is  unstable 
for  these  wave  lengths.  For  high  accuracy  machines  rounding 
errors  are  usually  very  small.  In  addition,  they  have  a 
random  character.  Accordingly,  one  will  obtain  acceptable 
solutions  even  though  short  waves  are  not  damped.  This 
optimistic  assessment  may  not  hold  for  nonlinear  problems. 

The  primary  reason  that  the  Courant  number  for  explicit 
finite  difference  methods  applied  to  hyperbolic  equations  must 
be  smaller  than  1  is  stability.  However,  this  restriction  is 
also  needed  from  the  point  of  view  of  accuracy,  because 
the  solutions  have  the  same  waviness  in  the  space  and  in  the 
time  directions.  If  one  admits  Courant  numbers  greater  than 
1,  then  one  loses  information  contained  in  the  approximation 
to  the  initial  conditions.  There  may,  however,  be  problems 
where  a  mesh  finer  than  required  for  reasons  of  accuracy  is 
used  in  the  space  direction. 

The  latent  presence  of  short  waves  makes  itself  felt 
in  the  semidiscretized  approach  during  the  integration  of  the 
resulting  system  of  ordinary  differential  equations.  In 
an  explicit  predictor  corrector  method  the  stability  limit 
is  about  the  reciprocal  value  of  the  largest  eigenvalue. 
Usually  this  limits  the  step  size  even  if  the  initial  con¬ 
ditions  do  not  contain  particular  solutions  pertaining  to  the 
large  eigenvalues  which  cause  instability  because  such 
particular  solutions  are  excited  by  the  truncation  errors  of 
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the  integration  scheme.  Implicit  methods  are  sometimes 
regarded  as  a  panacea,  but  this  is  not  entirely  correct. 
Basically,  one  deals  with  stiff  differential  equations. 
Familiar  solvers  for  stiff  differential  equations  are 
effective  for  problems  in  which  the  large  eigenvalues  which 
are  responsible  for  the  stiffness  have  large  negative  real 
parts  while  the  imaginary  part  is  fairly  small.  Under  these 
conditions  they  allow  the  use  of  large  integration  intervals 
in  regions  where  the  solution  is  sufficiently  smooth. 

Usually,  they  are  not  stable  if  the  real  part  of  the  critical 
eigenvalues  is  small  and  the  imaginary  part  is  large.  This 
is  the  present  situation. 

Our  discussions  have  been  restricted  to  rectangular 
elements  and  bilinear,  biquadratic  or  bicubic  shape  functions 
(This  makes  it  possible  to  discuss  space  and  time  dependence 
independently.)  The  accuracy  of  a  method  is  judged  by  the 
falsification  of  the  wave  speed  (and,  if  necessary,  also  of 
the  amplitude)  for  different  wave  lengths.  Results  of  this 
kind  are  shown  in  a  number  of  graphs.  If  one  uses  quadratic 
or  cubic  instead  of  linear  shape  functions,  then  one  has 
for  the  approximation  of  the  solution  in  the  space  direction 
twice  as  many  parameters  per  grid  point.  Each  discretization 
process  suppresses  waves  in  the  space  direction.  By  doubling 
the  number  of  parameters  per  grid  point  one  admits  wave 
lengths  down  to  one  half  of  the  original  limiting  wave  length 
In  this  sense,  a  representation  by  quadratic  or  cubic  shape 
function  is  equivalent  to  one  by  linear  shape  functions  with 
half  of  the  grid  size.  The  computational  effort  depends 
approximately  upon  the  overall  number  of  parameters.  We 
have  used  the  accuracy  of  the  wave  speeds  as  a  criterion  for 
a  comparison  of  different  methods.  The  mesh  size  for  linear 
shape  functions  is  chosen  one  half  of  that  for  quadratic  or 
cubic  shape  functions;  then  linear  and  cubic  elements  include 
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bicubic  shape  functions.  If  one  uses  bilinear  shape  functions 
and  weight  functions  of  the  same  kind  for  space  and  time, 
then  the  truncation  errors  cancel  if  the  Courant  number  is  1. 

The  method  is  unstable  if  the  Courant  number  exceeds  one. 

(An  exact  integration  of  the  differential  equations  for  the 
time  dependence  corresponds  to  Courant  number  zero.)  The 
favorable  results  obtained  with  Courant  number  one  have 
their  counterpart  in  finite  difference  methods. 

Occasionally,  particularly  in  the  development  of 
computational  procedures  for  the  transonic  problem,  it  is 
emphasized  that  one  must  use  a  difference  procedure  which 
reflects  the  marching  direction.  A  scheme  of  the  kind  just 
sketched  does  not  have  this  property;  a  marching  direction  is 
defined  by  the  way  in  which  the  initial  conditions  are 
prescribed . 

Seen  under  the  point  of  view  of  weighted  residual 
procedure  the  method  just  described  is  somewhat  unsatisfactory. 

In  computing  the  potential  at  a  new  time  station  it  considers 
not  only  the  residual  in  the  time  interval  that  is  newly 
computed,  but  also  in  the  preceding  interval  over  which  one 
no  longer  has  any  control . 

A  method  using  linear  shape  functions  has  been  investigated 
in  which  the  weight  function  is  constant  through  the  interval 
for  which  the  computation  is  carried  out  (say  from  t^  to  tfc+1 
with  an  infinitesimal  overlap  into  the  region  tk-1  to  tk-  The 
overlap  is  needed  in  order  to  take  into  account  the  delta 
function  in  the  second  derivatives  which  arises  at  the  point 
t^  because  of  a  jump  in  the  first  derivatives  which  is 
unavoidable  with  linear  shape  functions.  One  obtains  a 
method  which  is  stable  for  all  Courant  numbers,  the  solutions 
for  shorter  waves  are  rather  strongly  damped  and  one  is 
forced  into  rather  small  time  intervals  in  order  to  obtain 
an  acceptable  accuracy. 
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It  is  impossible  to  use  the  same  weight  functions  for 
the  space  and  the  time  directions  for  third  degree  shape 
functions  for  the  method  becomes  unstable  at  all  Courant 
numbers.  The  reason  is  rather  interesting.  For  third  degree 
shape  functions  one  obtains  the  solutions  at  time  t^+^;  that 
is,  the  values  of  and  u^+^  from  equations  which  contain 

the  (known)  values  of  u^,  u^,  u^_-^  and  u^_^.  In  order  to 
solve  the  differential  equation  in  the  interval  between 
tk  and  t^+1  it  suffices  if  one  known  u^  and  u^.  Speaking 
in  the  language  of  numerical  integration  techniques  for 
ordinary  differential  equations  one  would  say  that  the 
present  approach  includes  spurious  solutions.  These  are 
particular  solutions  of  the  homogeneous  discretized  system 
which  are  not  related  to  particular  solutions  of  the  original 
ordinary  differential  equation.  In  the  present  case  there 
are  two  such  spurious  solutions  and  one  of  them  is  always 
undamped  even  if  one  keeps  the  Courant  number  below  one. 

(For  the  other  particular  solutions  one  has  again  a  cancellation 
of  truncation  errors  for  space  and  time  if  the  Courant  number 
is  one.) 

One  can  devise  a  procedure  based  on  cubic  shape  functions 
in  which  no  spurious  solutions  are  encountered.  As  before, 
one  makes  the  transition  from  point  t^  to  point  t^+^  by  third 
degree  polynomials,  but  the  weight  functions  are  now  confined 
to  this  interval.  This  eliminates  the  objection  regarding 
the  use  of  the  residual  in  a  preceding  interval  to  determine 
the  solution  in  the  interval  under  investigation.  We  have 
studied  this  procedure  using  a  weight  function  1  over  an 
initial  part  of  the  interval,  namely,  t^  <  t  <  t^+c  ( ^+1  “  fck^  ' 

0  <  c  <  1,  and  again  over  the  remaining  part 
fck  +  c^tk+l  ”  fck^  <  <  ^+1*  Al1  particular  solutions  are 

unstable  for  c  <  1/2.  For  c  =  1/2  and  a  limited  Courant  number 
the  method  is  stable,  the  eigenfunctions  are  undamped.  There 
are,  however,  Courant  numbers  (somewhat  exceeding  one)  for 
which  unstable  particular  solutions  exist.  For  c  >  1/2  and 
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a  Courant  number  limited  approximately  to  a  value  below 
1  the  method  is  stable.  The  particular  solutions  are  slightly 
damped.  But  for  certain  higher  Courant  numbers  there  are 
always  unstable  particular  solutions.  The  procedure  is 
stable  for  all  Courant  numbers  in  the  limit  c  =  1.  Then 
the  weight  is  constant  throughout  the  interval  and  in 
addition  one  has  a  delta  function  weight  at  the  end  of  the 
interval.  This  amounts  to  a  mixture  of  finite  element  method 
and  collocation  method.  Even  in  this  limit  the  falsification 
of  long  waves  is  not  too  large. 

It  was  mentioned  above  that  from  the  point  of  accuracy 
Courant  numbers  larger  than  one  are  without  interest;  at 
least  if  the  shortest  wave  lengths  which  are  included  by  the 
discretization  in  space  direction  are  of  significance.  Cases 
where  one  will  admit  large  Courant  numbers  occur  in  transonic 
problems  because  they  include  the  vicinity  of  the  sonic  line. 
In  Section  VII  a  simple  problem  is  discussed  in  which  the 
Mach  number  assumes  the  value  one  at  one  boundary  of  the 
region.  In  the  direction  normal  to  the  flow  direction  an 
evenly  spaced  mesh  is  used.  For  each  length  of  the  interval 
one  obtains  a  maximum  value  of  y.  One  can  therefore  always 
choose  an  interval  in  the  downstream  direction  small  enough 
to  guarantee  stability.  If  the  mesh  size  in  the  direction 
normal  to  the  flow  is  multiplied  by  a  factor  a  and  one  wants 
to  maintain  the  same  Courant  number,  then  the  mesh  size  in 
the  downstream  direction  is  multiplied  by  a  factor  o'*.  In 
this  particular  problem  it  is  possible  to  find  particular 
solutions  by  a  product  hypothesis.  From  the  ordinary 
differential  equation  which  arises  for  the  direction  normal 
to  the  flow,  one  obtains  an  eigenvalue  problem  as  well  for 
the  ordinary  differential  equation  as  for  its  discretized 
form.  For  long  waves,  the  eigenfunctions  for  the  differential 
equation  and  for  its  discretized  form  are,  of  course,  very 
similar.  The  short  wave  eigenfunctions  are,  however,  entirely 
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unrelated.  For  the  differential  equation  they  are,  in 
essence,  given  by  Airy's  functions.  For  the  finite  difference 
formulation  the  highest  eigenfunctions  are  nearly  zero  (and 
oscillatory)  at  some  distance  from  the  boundary  which 
corresponds  here  to  Mach  number  1.  These  particular  solutions 
are  significant  only  at  the  first  few  mesh  points.  As 
contributions  to  the  solutions  of  the  original  partial 
differential  equation  problems  these  particular  solutions  are, 
of  course,  meaningless.  If  one  fails  to  observe  proper 
Courant  number  restrictions,  then  these  particular  solutions 
may  increase  as  one  moves  downstream.  The  falsification  so 
introduced  will  manifest  itself  in  some  roughness  in  the 
vicinity  of  the  line  corresponding  to  Mach  number  1,  but  at 
some  distance  from  it  the  error  will  be  very  small. 

In  the  examples  considered  up  to  this  point  the  weight 
functions  had  been  chosen  on  intuitive  grounds,  either 
because  of  their  resemblance  to  the  shape  functions  or,  for 
reasons  of  simplicity  as  constants.  Section  VIII  considers 
a  scalar  differential  equation  with  a  fixed  value  of  y  (in 
essence  the  reciprocal  of  the  wave  length) .  Then  one  can 
discuss  the  effect  of  a  residual  by  writing  down  the  solution 
of  the  inhomogeneous  differential  explicitly.  We  have  considered 
such  a  solution  in  the  form  of  the  Duhamel's  integral.  This 
can  be  regarded  as  a  solution  in  terms  of  the  Green's  function 
belonging  to  this  particular  one  dimensional  problem.  The 
integrals  occurring  here  have  the  form  of  weighted  residual 
expressions.  For  this  particular  problem  one  obtains  exact 
expressions  at  the  grid  points  if  the  solutions  of  the 
homogeneous  problems  are  taken  as  weight  functions.  In  reality, 
one  deals,  of  course,  simultaneously  with  a  whole  spectrum  of 
values  of  y.  At  best,  one  can  choose  weight  functions  which 
are  close  to  ideal  within  some  range  of  values  of  y.  But  the 
fact  is  of  interest  that  one  can  attune  the  weight  functions 
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to  a  certain  frequency.  In  an  example  we  have  chosen  y  =  0 
for  this  value.  The  weight  functions  are  then  given  for 
each  interval  by  a  constant  and  by  a  linear  function.  Prom 
our  numerical  results  we  cannot  claim  that  this  method  is 
superior  to  the  method  discussed  earlier,  but  the  author 
believes  that  the  theoretical  insight  is  valuable.  A  wider 
frequency  could  probably  be  covered  if  one  does  not  attune 
the  method  to  frequency  zero.  In  the  elliptic  region  one 
might  start  directly  from  the  Green ' s  function  and  impose 
the  requirement  that  the  lowest  order  terms  in  the  development 
of  the  Green's  functions  with  distance  vanish.  For  the 
Laplace  equation  this  gives  linear  weight  functions.  For  the 
Helmholtz  equation  there  are  always  some  contributions  of 
the  lowest  order  terms  which  vanish  with  distance  only  weakly. 
They  become  small  only  if  the  mesh  is  sufficiently  fine. 

The  same  idea  applied  to  hyperbolic  problems  fails, 
for  the  hyperbolic  distance  is  zero  for  points  that  lie  on 
the  characteristics  through  the  points  where  the  residuals 
occur.  The  weight  functions  for  hyperbolic  problems  used 
in  the  preceding  analysis  are  suitable  for  long  waves  in  time 
and  space.  Short  waves  are  treated  inaccurately,  at  best. 

A  finite  element  of  this  kind  is  therefore  unsuited  to  treat 
discontinuities  which  propagate  along  characteristics. 

For  two  dimensional  problems  it  is  possible  to  combine 
the  idea  of  finite  elements  with  the  concept  of  characteristics. 
This  is  done  in  Section  IX.  The  shape  functions  are  then 
defined  in  characteristic  quadrangles.  Important  is  the 
result  that  the  weight  functions  should  be  constant  in  one 
of  the  characteristic  directions.  With  this  choice  it 
becomes  possible  that  contributions  of  residuals  which  lie 
on  the  same  characteristic  cancel  each  other.  The  desirability 
of  this  form  of  weight  functions  becomes  obvious  again  if  one 
considers  the  Green's  function  of  this  problem.  A  finite 
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element  method  carried  out  in  this  manner  differs  only 
slightly  from  the  method  of  characteristics.  Within  an 
element  discontinuities  are  suppressed,  but  along  the  element 
boundary  discontinuities  are  admitted  because  they  do  not 
contribute  to  the  residual. 

The  method  of  characteristics  is  not  easily  carried 
over  to  three  dimensions  and  no  conclusions  in  this 
direction  have  been  reached. 

The  picture  that  emerges  from  these  discussions  is 
somewhat  discouraging.  A  Courant  number  limitation  always 
exists,  if  not  for  reasons  of  stability  then  for  regions 
of  accuracy.  Higher  (third)  order  elements  are  advantageous 
because  they  give  for  long  wave  lengths  a  much  better 
representation  of  the  wave  speed  than  linear  elements. 

The  use  of  third  order  splines  is  probably  advantageous 
because  it  suppresses  short  waves  which  have  an  undesirable 
effect  on  the  time  integration  and  still  give  a  good 
representation  for  long  waves.  The  method  of  characteristics 
possibly  in  combination  with  the  finite  element  concept  may 
still  prove  a  superior  alternative. 

The  report  goes  somewhat  further  than  an  evaluation  of 
different  methods  by  means  of  sample  problems,  it  relates 
the  results  that  could  be  found  in  this  manner  to  the 
behavior  of  particular  solutions  of  the  approximating 
equations.  In  this  manner  certain  inherent  difficulties 
directly  connected  with  the  nature  of  hyperbolic  equations 
become  obvious.  The  approaches  discussed  here  may  not  be 
the  last  word  in  this  problem  area.  Other  more  ingenious 
approaches  may  well  exist,  but  even  they  will  encounter  the 
same  obstacles  and  then  a  familiarity  with  these  phenomena 
may  be  valuable. 
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APPENDIX  I 

THE  VARIATIONAL  FORMULATION  OF  REDDY 


A  variational  formulation  for  some  hyperbolic  problems 
which  include  the  boundary  conditions  in  a  proper  manner  has 
been  given  by  Reddy  (Ref.  1) .  The  author  believes  that  a 
variational  principle  by  itself  (in  contrast  to  an  extremum 
principle)  is  not  particularly  useful  for  numerical  work. 

The  present  appendix  shows  in  a  simple  example  the  essence  of 
Reddy's  formulation.  Consider  the  one-dimensional  problem 

*»'•*/.**■  -  fft‘-  »  (A 


with  initial  conditions  of  the  type  encountered  in  a  hyperbolic 
problem 

u,{0)  -  ^  J  to(O)  -  u!0  (A2) 

The  independent  variable  is  t,  the  derivative  of  a  function  with 
respect  to  its  argument  is  denoted  by  prime.  Let  the  final 
value  of  t  be  tQ  >  0. 

Now  consider  the  functional 

ytvc)  •  -  uJ (0)n,(tc  )  +  f[i  uZ(?)U(to-t)  ?0  -V )  (A3 } 

*  -  /(r)to  fo 

We  form  the  variation  of  i|>  (u)  ,  postulating  that  the  functions 
u  in  competition  satisfy  eqs.  (A2) .  Then  one  has 


$  io(o)  •  o  .  «P  ul(o)-  o 


(A4) 


One  obtains  to 

Syfu,)  m  -lofo)  £u,(tj  +f {if *(Vlu!fo-Z)  -tj 

+if0  -VJ+1& 

-  Stofa-trjJi/tr 


One  obtains  by  carrying  out  an  integration  by  parts  in  the  first 
two  terms  of  the  integrand 

fy/tc)  .  -to/c)ftofc)  +£[ rtuft)tcfa-v)-  <vft) }*(■%- 

+ J'fl'fM'ftJlSft -t)  +i  toft)  tufa  -t) 

*  i  ft  tic(V)Utft*  +if»  uft)  tufa- 1)  -fft)tufa 

In  the  first  and  third  terms  of  the  integrand  the  umbral  variable 
t  is  replaced  by  tQ  -  x.  In  addition,  we  evaluate  the  terms 
outside  of  the  integral 

tf/u)  *  -u!fo)tufa)  Siofajtofa)  -  tic/ejuffa) 

/ 

The  terms  outside  of  the  integral  vanish  because  of  the  boundary 
conditions  Eqs.  (A2)  and  (A4) .  By  the  requirement  that  the 
variation  of  6ijj(u)  vanish  one  obtains  indeed  Eq.  (Al)  . 

The  last  expression  can  be  interpreted  as  a  weighted 
residual  formulation.  The  weight  functions  are  n(t)  =  6u(tQ-t). 
Let  the  shape  functions  be  restricted  to  some  subspace  of  the 
space  of  admissible  functions.  Strictly  speaking  one  cannot 
maintain  that  the  shape  functions  n(t)  belong  to  the  same 
function  space,  for  the  shape  functions  and  their  first 
derivatives  vanish  for  argument  0,  the  shape  functions  and 
their  derivative  vanish  for  argument  t q.  The  property 
characteristic  of  the  usual  variational  formulations,  namely 
that  weight  functions  and  shape  functions  are  taken  from  the 
same  subspace  no  longer  applies.  The  method  would  probably 
encounter  difficulties  if  one  uses  an  uneven  grid  in  t.  One 
notices  furthermore  that  the  procedure  requires  that 
p(t)  =  p{tQ-t)  (in  our  example  we  have  assumed  pQ  =  const). 

These  are  obstacles  to  the  application  of  Reddy's  formulation 
to  numerical  work. 


J  [toft)  *  ft  -/W  fa -VS* 
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APPENDIX  II 


OBSERVATIONS  REGARDING  THE  TREATMENT  OF 
BOUNDARY  CONDITIONS  IN  THE  METHOD  OF  WEIGHTED  RESIDUALS 


The  observations  made  in  this  appendix  arose  from  dis¬ 
cussions  which  I  had  with  Dr.  Soliman  about  the  method  of 
weighted  residuals.  Although  they  do  not  refer  specifically 
to  hyperbolic  problems  they  are  included  here  because  they  serve 
to  round  out  the  knowledge  about  the  method  of  weighted 
residuals.  The  essential  idea  is  already  seen  in  a  one¬ 
dimensional  problem.  Consider  the  differential  equation 

-j%)  +  +/(*)•  0  (A6) 

with  boundary  conditions 


<fr(o)  *  0 

C  <j>(L)  *  ^(L)  *  $  -  0 


(A7) 

(A8) 


Including  the  second  boundary  condition  in  a  weighted  residual 
formulation  we  write 

L  *  f 

J [-<}(*)  t  fftxjft)  *  f(i)]f(xu/t  +  ( [c,  %]*  o  (A9 

9 

where  4)  is  a  Weight  function  and  C  is  a  single  weight.  The 
functions  $  in  competition  are  assumed  to  satisfy  the  boundary 
condition 


We  postulate 

ffO)*  0 

The  requirement  of  existence  of  second  derivatives  in 
<)>  can  be  relinquished  if  one  carries  out  an  integration  by 
parts.  One  obtains,  after  substitution  of  the  boundary  conditions 

ft y<k) /ft)  +fMfft) fM  *  ftoffKpx  -  V  c4) 

m  O  (A10) 
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If  one  chooses  C  independent  of  the  weight  function  tj>,  then 
one  obtains  directly  the  condition  Eq.  (A8) ;  the  admissible 
functions  4>  must  satisfy  the  boundary  conditions  at  x  =  L  as 
well  as  x  =  0. 

According  to  Soliman,  one  can  simplify  the  problem  by 
choosing 

^Ct.  *  (All) 

with  this  choice  4>'(L)  vanished  in  Eq.  (A10)  .  This  may  well 
be  worthwhile  particularly  in  multidimensional  problems 
because  it  obviates  the  need  to  evaluate  the  derivative  in 
the  direction  of  the  conormal  to  the  boundary. 

To  justify  this  special  choice  we  rewrite  Eq.  (A9) 
expressing  C  by  Eq.  (All) .  Then  one  has 

f {-<f%) +?(*)/(*) +/(*)+&*  +</>'(*>+  irIjwJfO 

0  0 

*>0  (A12) 

Here  the  term  generated  by  the  boundary  condition  is  included 
in  the  integral  by  means  of  the  delta  function  6(x,  L-e)  . 

The  6  singularity  approaches  the  point  x  =  L  from  the  inside 
of  the  interval  as  e  +  0.  One  sees  that  the  failure  to  satisfy 
the  boundary  condition  at  x  =  L  appears  in  the  weighted 
integral  for  the  residual  in  the  form  of  a  delta  function. 

There  are  other  points  within  the  interval  where  the  residual 
is  allowed  to  have  delta  singularities,  these  are  points  where 
<p '  has  jumps.  The  special  choice  of  C  in  Eq.  (All)  therefore 
has  the  effect  of  balancing  the  failure  to  satisfy  the  boundary 
conditions  against  other  residuals  caused  in  the  differential 
equation  by  the  approximation  to  <(>.  This  is,  of  course  legitimate. 
The  situation  is  similar  for  multidimensional  problems.  The 
boundary  terms  in  Eq.  (A12)  are  then  replaced  by  integrals 
over  the  boundary  surface. 
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The  factor  by  which  the  delta  function  is  multiplied 
tends  to  infinity  if  c2  tends  to  zero.  Then  this  approach 
must  be  abandoned.  In  this  case  the  functions  $  must 
satisfy  the  boundary  conditions  at  x  -  L  exactly. 


In  the  problem  at  hand  one  has  also  a  variational 
formulation ,  namely 

tf/f-rSrt*  ”  r  £#11*  &(L,}m  0 

In  the  case  c2  -  0 ,  that  is,  if  the  value  of  <t>  is 
prescribed  at  the  boundary,  one  can  proceed  in  different 
ways.  In  one  formulation  one  postulates 
L  ,  t 

S  f  [{  i  *  0 


(A13) 


and  postulates  that  the  admissible  functions  $  satisfy  the 
boundary  conditions  at  x  =  0  and  x  =  L.  Alternatively,  one 
can  introduce  the  boundary  conditions  at  x  =  L  by  means  of  a 
Lagrange  multiplier  and  disregard  the  boundary  condition  at 
x  =  L  during  the  variations.  This  formulation  is  discussed 
because  it  has  some  similarity  to  Eq.  (A10) .  However,  we  shall 
see  that  the  two  formulations  are  different  in  principle. 


One  has 


£[ J[i </!(*) +ieftyftx> * f &]*  0 


(A14) 


Hence 


/  /y  4  ffy))tyt*)dx  4-  (4'fa  +VfyfL)  ■  0 

« 

and 

-4%  +{(*>•  0  (A15) 

^  *(L) 

The  variation  vanishes  if  one  solves  the  original  differential 
equation  with  the  boundary  conditions  <M0)  =  0  and  $'(L)  «  -  X. 


This  can  be  done  for  any  value  of  X.  The  desired  solution 
is  obtained  if  one  determines  X  so  that  the  original  boundary 
condition 

A//  )  _  _  A 


is  satisfied. 


In  a  weighted  residual  approach  one  would  write 

f  f -$*(*)  < f(U  fa)*  ftl>)  pfV  dx  •  o 

0 

$/c,)  *  0 


(A16) 


Eqs.  (A16)  and  (A15)  are  not  identical.  In  the  approach  with 
Langrange  multipliers  one  computed  in  principle  a  family  of 
solutions  which  satisfy  a  different  boundary  condition  specified 
by  the  choice  of  X.  From  this  family  one  then  selects  the 
one  which  satisfies  the  boundary  condition  actually  prescribed. 
In  a  weighted  residual  approach,  the  boundary  conditions  are 
introduced  directly. 

Functions  <J>  which  fail  to  satisfy  the  boundary  conditions 
give  approximations  with  a  jump  of  <J>  at  the  boundary.  Such  a 
jump  could  be  accommodated  by  means  of  the  derivative  of  a 
6  function.  But  such  derivatives  cannot  be  dealt  with  by  a 
method  of  weighted  residuals.  This  is  important  in  multi¬ 
dimensional  problems  where  it  is  frequently  impossible  to 
find  shape  functions  <j)  which  satisfy  the  boundary  conditions 
exactly.  It  is  preferable  to  consider  the  error  introduced 
by  the  failure  to  satisfy  the  boundary  conditions,  by  itself, 
rather  than  to  attempt  to  include  it  in  an  overall  weighted 
residual  formulation. 
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APPENDIX  III 

CHECK  EXPRESSIONS  FOR  SOME  FORMULAE  OF  SECTION  II 


The  Eqs.  (11)  and  (13)  of  cases  2a  and  2b  respectively 
can  be  checked  by  substituting  the  values  of  ,  <J>k  and 

4>k+^  for  three  piecewise  linear  functions.  For  such 
functions  the  operator  on  the  left  of  these  equations  must 
give  the  expressions  which  one  obtains  by  substituting  them 
and  the  chosen  weight  functions  into  Eq.  (7)  . 

The  following  expressions  will  be  chosen 


Hi,t)  .  Cft) 


Hi.  v 


UcH) 

h 

f  O 


to  t  ay  <0 

'f)£  )  *  £  ay/h  Cft)  &y>0 


(Al7a) 


(A17b) 


(A17c) 


In  Case  2a  the  weight  functions  are  given  by  Eq.  (10)  .  One 
obtains  by  substituting  Eqs.  (10)  and  (Al7a)  into  Eq.  (7) 


J  f*y ft)  ~  /  0  Hf&ylbidfajl1*)  »  (A18) 

The  expression  (Al7b)  substituted  into  Eq.  (7)  gives  zero 
2  2 

because  3  <b  and  3  are  antisymmetric  and  w(Ay)  is  symmetric 

3?.  3y 

The  expression  (Al7c)  substituted  into  the  first  term 
of  Eq.  (7)  gives 

£ /(*/*) 

o 

The  second  term  must  be  evaluated  in  the  sense  of  generalized 
functions  /  .. 
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■ 


One  obtains,  therefore,  for  case  (2b)  with  the  "test"  function 
Eq.  (A17c) 

+  (A19, 

One  finds  from  (Al7c) 


k-t  ~ 


A 


4> 


A.,  * c(t  ‘ 


(A20) 

One  obtains  the  expression  (A19)  if  one  substitutes  Eq.  (A20) 
into  Eq.  (11) . 


Case  2b. 

With  the  weight  given  by  Eq.  (12)  one  obtains  in  the 
same  manner 

d2C/dt2  for  Eqs .  (A. 17a)  and  (12) 

0  for-Eqs.  (A. 17b)  and  (12) 

d2C/dt2  -  h  2C  for  Eqs.  (A. 17c)  and  (12) 


Cases  3. 

In  the  resulting  formulae,  Eqs.  (16)  ,  (18)  and  (20) 
there  appears  six  unknowns,  ^^(t),  (^^(t),  <fk(t),  4>£(t), 
<J>k+1(t)  and  For  a  check  of  these  equation  one  needs 

six  independent  functions  <f>,  which  are  piecewise  of  third 
degree  and  continuous  in  function  and  first  derivative.  One 
notices  that  the  first  and  second  equations  of  the  equation 
pairs  (16)  ,  (18)  and  (20)  are  automatically  satisfied  if  <j> 
is  antisymmetric  and  symmetric,  respectively  with  respect  to 
the  point  y  =  y^. 

A  convenient  choice  of  the  test  functions  is 


<j>  =  C(t)  (A21a) 

4>  =  C  ( t)  Ay/h  (A21b) 

4>  -  C(t)  (Ay/h)2  (A21c) 

<j>  -  C(t)  (Ay/h)3  (A21d) 


96 


(A21e) 


f-C(t) (Ay/h)3  Ay/h  <  0 

L  C(t)(Ay/h)J  Ay/h  >  0 

(-C(t)  (Ay/h)  2  Ay/h  <  0 

♦«  j  ,  (A21f) 

L  C(t) ( Ay/h) z  Ay/h  >  0 


Case  3a. 

Substituting  the  expression  (A21)  together  with  the  weight 
functions  (15)  into  the  operator  on  the  left  of  Eq.  (7) ,  one 
obtains 


dzC/dtz 

for 

Eqs . 

(A21a) 

and 

(15a) 

0 

for 

Eqs . 

(A21a) 

and 

(15b) 

0 

for 

Eqs. 

(A21b) 

and 

(15a) 

1/15  d2C/dt2 

for 

Eqs. 

(A21b) 

and 

(15b) 

2/15  d2/dt2  -  2h~2C 

for 

Eqs. 

(A21c) 

and 

(15a) 

0 

for 

Eqs. 

(A21c) 

and 

(15b) 

0 

for 

Eqs . 

(A21d) 

and 

(15a) 

2  d2C  -2.-2 
105^1  5 " 

for 

Eqs. 

(A21d) 

and 

(15b) 

i  £h-2c 

dt  3 

for 

Eqs. 

(A21e) 

and 

(15a) 

0 

for 

Eqs. 

(A21e) 

and 

(15b) 

0 

for 

Eqs . 

(A21f) 

and 

(15a) 

J_  -  ih-2c 

70  7“  U 

for 

Eqs. 

(A21f) 

and 

(15b) 

ar 


These  results  can  be  used  to  check  Eqs.  (16) 
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Case  3b. 

The  weight  functions  axe  given  by  Eqs.  (17),  the  test 
functions  by  Eqs.  (A21) .  The  following  results  are  obtained 
by  substitution  into  Eq.  (7) 


2  2 
aC/dz 

for 

Eqs . 

(A21a) 

and 

(17a) 

0 

for 

Eqs . 

(A21a) 

and 

(17b) 

0 

for 

Eqs. 

(A21b) 

and 

(17a) 

\  d2C/dt2 

for 

Eqs. 

(A21b) 

and 

(17b) 

Y2  d2C/dt2 

-  2h_2C 

for 

Eqs. 

(A21c) 

and 

(17a) 

0 

for 

Eqs . 

(A21c) 

and 

(17b) 

0 

for 

Eqs. 

(A21d) 

and 

(17a) 

h  d2c/dt2 

-  2  h  C 

for 

Eqs. 

(A2  Id) 

and 

(17b) 

22  d2C/dt2 

-  |  h“2C 

for 

Eqs. 

(A21e) 

and 

(17a) 

0 

for 

Eqs. 

(A21e) 

and 

(17b) 

0 

for 

Eqs. 

(A21f) 

and 

(17a) 

j2  d2C/dt2 

-  2h_2C 

for 

Eqs. 

(A21f ) 

and 

(17b) 

These  results  can  be  used  to  check  Eqs.  (18) . 

Case  3c. 

The  weight  functions  are  given  by  Eqs.  (19) .  The  test 
functions  are  given  by  Eq.  (A21) .  The  following  results  are 
obtained  by  substitution  into  Eq.  7. 

j  d2C/dt2  for  Eqs.  (A21a)  and  (19a) 

2  d2C/dt2  for  Eqs.  (A21a)  and  (19b) 
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0 

for 

Eqs. 

(A21b) 

and 

(19a) 

I  d2C/dt2 

for 

Eqs. 

(A21b) 

and 

(19b) 

h  d2c/dt2 

-  h"2C 

for 

Eqs. 

(A21c) 

and 

(19a) 

13  A2n/A  +  2 

96  a  C/dt 

-  h“2C 

for 

Eqs. 

(A21c) 

and 

(19b) 

0 

for 

Eqs. 

(A21d) 

and 

(19a) 

it  a2c/dt2 

- 1  h*2c 

for 

Eqs. 

(A21d) 

and 

(19b) 

512  d2c/dt2 

-  (3/8) h”2C 

for 

Eqs. 

(A21e) 

and 

(19a) 

is  d2c/dt2 

-  §  h"2c 

for 

Eqs. 

(A21e) 

and 

(19b) 

0 

for 

Eqs. 

(A2  If) 

and 

(19a) 

13  A2n/A+2 
jjg-  d  C/dt 

-  h“2C 

for 

Eqs. 

(A21f) 

and 

(19b) 

results  are 

used  to  check  Eq. 

20. 

Cases  4. 

The  resulting  formulae,  Eqs.  (21)  and  (23)  are  correct 
for  continuous  piecewise  quadratic  functions.  In  Eqs.  (22) 
and  (24)  there  appear  five  different  functions.  Therefore, 
one  needs  five  test  functions.  The  following  expressions 
are  suitable. 


4>  =  c(t) 

(A22a) 

<J)  =  C(t)  Ay/h 

(A22b) 

4>  =  C(t)  (Ay/h)2 

(A22c) 

r-C(t)  Ay/h 

*  - 

^  C(t)  Ay/h 

-1  <  Ay/h  <  0 

0  <  Ay/h  <  1 

(A22d) 
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♦  - 


-  -C(t)  (Ay/h)2 

■ 

.  C (t) (Ay/h) 2 


-1  <  Ay/h  <  l 
0  <  Ay/h  <  1 


(A22e) 


There  is  no  need  to  examine  the  test  function  Eq.  (A22e) 
separately.  Weight  functions  (21a)  and  (23a)  gives  zero 
for  reasons  of  symmetry,  for  weight  functions  Eqs.  (21b) 
and  (23b)  it  coincides  with  case  (A22c)  because  the  weights 
are  identically  equal  to  zero  outside  of  the  region 
0  <  Ay/h  <  1.  The  following  results  are  obtained  by 
substitution  into  Eq.  (7) . 


Case  4a. 

The  weight  functions  are  given  by  Eqs.  (21)  the  test 
functions  by  Eqs.  (A22) .  The  following  results  are  obtained 
by  substitution  into  Eq.  (7) . 


d2C/dt2 

for  Eqs. 

(A22a) 

and 

(21a) 

|  d2C/dt2 

for  Eqs. 

(A22a) 

and 

(21b) 

0 

for  Eqs. 

(A22b) 

and 

(21a) 

12  2 
|  aC/dz 

for  Eqs. 

(A22b) 

and 

(21b) 

^  d2C/dt2  -  2h_1C 

for  Eqs. 

(A22c) 

and 

(21a) 

^  d2C/dt2  -  j  h"2C 

for  Eqs. 

(A22c) 

and 

(21b) 

j  d2C/dt2  -  2C 

for  Eqs. 

(A22d) 

and 

(21a) 

j  d2C/dt2 

for  Eqs. 

(A22d) 

and 

(21b) 

The  results  are  used  to  check  Eq.  (22). 
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Case  4b 


The  weight  functions  are  given  by  Eqs.  (23)  and  the 
test  functions  by  Eqs.  (A22) .  The  following  results  are 
obtained  by  substitution  into  Eq.  (7). 


\  d2C/dt2 

for  Eqs. 

(A22a) 

and 

(23a) 

\  d2C/dt2 

for  Eqs. 

(A22a) 

and 

(23b) 

0 

for  Eqs. 

(A22b) 

and 

(23a) 

i  d2C/dt2 

for  Eqs. 

(A22b) 

and 

(23b) 

^  d2C/dt2  -  h‘2C 

for  Eqs. 

(A22c) 

and 

(23a) 

||  d2C/dt2  -  h"2C 

for  Eqs. 

(A22c) 

and 

(23b) 

d2C/dt2  -  2h~2C 

for  Eqs. 

(A22d) 

and 

(23a) 

|  d2C/dt2 

for  Eqs. 

(A22d) 

and 

(23b) 

The  results  are  used  to  check  Eq.  (24) . 
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APPENDIX  IV 


DIAGONAL  FORM  OF  A  SYSTEM  OF 
ORDINARY  DIFFERENTIAL  EQUATIONS 


The  system  of  differential  equations  whose  stability 
is  studied  is  given  by 

+  M2  +  *  0  (A23) 

-+■  '*-*■'* 

Here  y  and  f  denote  respectively  the  dependent  variable 
and  a  known  vector  with  n  components.  L  and  M  are  n  by 
n  matrices.  It  is  assumed  that  the  determinant  of  L  does 
not  vanish.  In  writing  our  equations  we  consider  the 
vectors  as  n  by  1  or  1  by  n  matrices.  The  dot  denotes 
differentiation  with  respect  to  the  independent  variable  t. 

The  matrices  L  and  M  are  independent  of  t.  (This  assumption 
is  always  made  in  stability  discussions.  Heuristically ,  it 
is  justified  by  the  fact  that  for  an  arbitrarily  fine  mesh 
in  an  equation  with  variable  coefficients  one  can  carry  out 
a  considerable  number  of  integration  steps  before  the  matrices 
L  and  M  change  appreciably.)  We  associate  with  Eq.  (23)  the 
generalized  eigenvalue  problem 

(H  ~  s  0  (A24) 

and  its  adjoint 

'  °  (A25) 

where  Tk  is  a  n  by  1  matrix  (i.e.  a  column  vector)  and  is 
a  1  by  n  matrix  (a  row  vector),  X.  and  X.  are  respectively 

til  til  ^  ^ 

the  k  and  £  eigenvalue.  Proceeding  in  a  familiar  manner 
one  has 


Hence 


•i  t 
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Let  T  and  S  be  n  by  n  matrices;  the  k  column  of  T  is  given 
by  T^,  the  row  of  S  by  S^.  Assume  that  S  and  T  are 

normalized  so  that 


(A26) 

where  'I  is  the  n  dimensional  identity  matrix.  Then  one 
obtains  from  Eq.  (A24) 


StlT  -  A 


(A27) 


where  A  is  a  diagonal  matri  :  whose  kth  element  is  Now 

set 


Ss  7'u 

y -  tv 

where  u  is  a  column  vector  with  n  components.  Then  one 
obtains  from  Eq.  (A23)  and  by  premultiplication  with  S  and 
by  using  Eqs.  (A26)  and  (A27) 


(A28) 


tO  *  O  (A29) 

This  is  a  system  of  single  equations,  for  4  is  a  diagonal 
matrix. 


This  fact  allows  one  to  discuss  the  stability  of  the 
system  in  terms  of  the  stability  of  a  single  scalar  equation. 
Specifically,  we  shall  derive  the  fact,  familiar  from  the 
theory  of  integration  process,  that  one  obtains  the  same 
result  if  one  first  diagonalizes  the  system  of  equations  and 
then  applies  a  predictor  corrector  method  or  first  applies 
the  predictor  corrector  method  and  then  diagonalizes.  The 
same  holds  for  a  finite  element  procedure  applied  to  the 
time. 


In  a  predictor  corrector  method,  no  matter  which  specific 
scheme  one  applies,  one  always  has  as  predictor  formula 
fo) 

-Ji  i  *  ••  • 


j 


103 


The  values  of  the  vector  y  after  the  l  corrector  iteration 
is  characterized  by  a  superscript.  The  values  of  y  and  y*  at 
the  point  i  and  preceding  points  are  fixed,  therefore  no 
reference  to  an  iteration  step  is  needed.  One  has  in  the 
corrector  step 

&  >m  &  *  *  £  &  *  A  yt-2. y 

The  a^'s,  0j's,  6j's  and  d^'s  are  constants  independent  of 
i.  In  addition,  one  has  Eq.  (A23)  to  be  satisfied  at  the 
individual  grid  points.  Substituting  Eq.  (A28)  into  the 
last  two  equations  and  premultiplying  by  T-1,  one  obtains 
the  same  expressions  but  with  y  and  y‘  replaced  by  u  and  u*. 

In  addition,  one  has  Eq.  (A29) .  In  the  resulting  equations 
the  vectors  u^  and  uV  are  multiplied  by  constants  or  by  the 
diagonal  matrix  A.  The  system  can  therefore  be  decomposed 
into  its  individual  components.  The  further  discussion  can 
therefore  be  carried  for  each  component  in  Eq.  (A29)  separately 
each  with  only  one  independent  variable.  The  argument  for 
implicit  methods  is  the  same. 

To  apply  a  finite  element  method  to  the  solution  of 
Eq.  (A23)  one  multiplies  this  equation  with  scalar  weight 
functions  wm(t)  (usually  of  finite  support)  and  integrates 
with  respect  to  t. 

ft. 

+  f  H*, fr jJt  +  o 

1/  t/  ^ 

Now  we  substitute  Eqs.  (A28) ,  premultiply  the  resulting 
equation  by  S  and  apply  Eqs.  (A26)  and  (A27) .  One  obtains 


+  f  o  (A3Q) 

-t,  i,  *• 

This  expression  is  identical  with  the  one  which  one  would 
have  obtained  by  applying  the  weight  function  wm(t)  directly 
to  Eq.  (A29) .  Since  A  is  a  diagonal  matrix,  one  can  discuss 
each  component  of  u  separately.  This  result  is  the  basis  of 
the  analysis  in  Section  VI. 
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APPENDIX  V 

FUNDAMENTAL  SOLUTIONS  AND  RELATED  SUBJECTS 


In  this  appendix  the  wave  equation  is  regarded  as 
the  linearized  equation  for  a  supersonic  flow  at  Mach 
number  /l.  We  consider  accordingly 


(A31) 


The  perturbation  velocities  are  then  given  by  the  vector 
grad  <}>.  We  determine,  also,  the  perturbation  of  the  mass 
flow  vector.  In  the  nonlinearized  problem  the  mass  flow 
vector  is  given  by  pgrad$,  where  p  depends  upon  grad  <|> . 


Specifically,  one  has  from  Bernoullis,  equation 
*  ?  Jfi  {&+<&})  -  O 

and  therefore,  with  the  sound  velocity  "a"  given  by 

</f*  -  +4?^) 

Now  the  total  potential  4  is  given  by 

»  Ux  *  $ 

Hence,  for  the  x  component  in  the  perturbation  of  the  mass 
flow 

f'fxfi-  u/a-) 

and  for  the  y  component 

Here  pQ  is  the  density  in  the  unperturbed  flow.  One  therefore 
obtains  for  the  x  and  y  components  of  the  perturbation  of 
the  mass  flow  vector  at  Mach  number  S2,  respectively  -P04>x 
and  pQ<J>y*  (Important  is  the  negative  sign  in  the  x  component.) 
The  potential  equation  in  the  form 

Ax  *  jytj  )  ■  0 

expresses  the  conservation  of  mass.  By  the  usual  integration 
by  parts  (Gauss'  theorem)  one  obtains 

*  "A  (A32) 
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where  one  travels  along  the  path  of  integration  in  the 
direction  for  which  the  region  R  is  to  the  left. 

Assume  that  some  portion  of  the  contour  is  given  by  a 
left  going  characteristic  y  =  c  +  x  .  Then  the  correction 
to  the  mass  flow  through  this  part  of  the  contour  is  given 

by  -A  Jff,  * 


where  ds  is  the  line  element.  In  this  case 

Thus,  one  obtains  as  a  correction  to  the  mass  flow  passing 
through  a  left  going  characteristic 

"  A  /  7^  *  -  A  fa  ) 

Analogously,  for  a  right  going  characteristic 


(A33) 


-  J>e  fffi  T-dl  m  hffoj'fori 


(A34) 


Eq.  (A31)  is  solved  by 

.  /,(y+x)  + 

To  obtain  a  fundamental  solution  we  postulate  that  $  =  0  for 
x  <  0  and  that  at  x  =  0,  <|>x  is  given  by  the  negative  of  a 
delta  function  at  the  origin.  The  requirement  <f>  s  0  at  x  =  0 
is  satisfied  by  choosing 

f  *  ~  f  =  f 

Then  one  has 

&  m  ffy+x*  t  ffy~x) 

It  follows  that 

f(*)a  O  l<<> 

f (*)»'£ 

Thus,  one  finds  for  x  >  0  as  the  potential  caused  by  a  source 
at  the  origin 

4>  =  0  y  <  -x 

<j>  -  -1/2  -x  <  y  <  x  x  >  0  (A3, 

<p  =  0  y  >  x 

« 

<p  -  0  x  <  0 


(A35) 
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If  one  crosses  the  right  going  characteristic  y  =  -x  along 
a  left  going  characteristic,  then  the  correction  to  the 
mass  flow  is  according  to  Eq.  (A33)  (1/2) Pq.  The  same  result 
is  obtained  if  one  corrects  the  characteristic  y  =  x.  A 
fundamental  solution  gives  a  correction  to  the  mass  flow 
vector  along  the  characteristics  through  the  origin  in  the 
direction  of  the  characteristics  in  the  form  of  a  delta 
function  of  intensity  l/2pQ.  Everywhere  else  the  correction 
is  zero  (Figure  30.) .  The  gradient  of  the  velocities  is 
perpendicular  to  these  characteristics.  The  perturbation 
does  not  die  out  with  distance. 

Except  for  a  constant,  the  fundamental  solution  for 

2  2 

the  Laplace  equation  is  given  by  log(x  +  y  ) .  The  corresponding 
expression  for  the  wave  equation  is  given  by  log(x  -  y  ). 

For  x  /  0,  y  /  0  and  x  ^  jyj  this  expression  certainly  satisfies 
the  partial  differential  equation,  but  obviously  it  does  not 
give  the  desired  fundamental  solution. 

Also  of  interest  is  the  perturbation  caused  in  an 
axisymmetric  flow  by  a  source  at  the  origin.  It  provides  the 
fundamental  solution  for  the  three  dimensional  wave  equation. 

The  linearized  equation  for  three  dimensional  supersonic  flows 
for  Mach  number  /2*  is  given  by 

*  fai  *  T  (A3 6) 

where  f(x,y,z)  is  the  local  source  strength.  One  obtains, 
by  integration  over  a  volume  R 

Jiff-fax  +faf  tjuMyS* 

The  right  hand  sides  gives  the  combined  strength  of  the 

sources  within  the  volume  R.  Now  we  assume  axial  symmetry 

and  introduce  cylindrical  coordinates 

x  -*  x 

r*-*y*+** 

One  then  obtains 

x  4 


(A37) 


i 


It  is  convenient  to  introduce  coordinates  oriented  according 
to  the  directions  of  the  characteristics  (Fig.  31) 


f-r+* 


V 


-r+x 


izJL 


*  t 

with  the  restriction 


The  x  axis  is  given  by  £  =  n • 

One  obtains  for  the  Jacobian  of  the  transformation 


d(g,n) 

d(r,x) 


2 


Moreover 

27  m  7J~  2 1  >2x^1  W 

Thus,  one  obtains  the  following  expression  for  the  total 

mass  flow. 

J. <A38’ 

In  these  coordinates  the  original  equation  for  the 
perturbation  potential ,  namely 
_ 

assumes  the  form 

&  *  klt!'v  *  (A39) 

Here  one  has  the  particular  solution 

•>/*.  /  */,k 

<f>  -  (X*-  rLJ  *1  (A40) 

This  is  easily  verified. 

The  expression  Eq.  (A38)  can  be  transformed  into  a  contour 
integral 

The  integrand  vanishes  at  the  x  axis  (5  =  ri)  • 
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To  determine  the  total  source  strength  we  choose  for 
R  a  region  bounded  by  the  x  axis  (£  =  n)  *  by  a  line  Z  s  Zi 
and  closed  to  the  left  by  some  arbitrary  curve.  Then  one 
finds 

i A/,  -i)  Ms?  - +  *Jr 

Substituting  the  specific  expression  for  <j>,  Eg.  (A40)  one 


obtains  .  .  i  -'A  V  ,  -AA  . 

**!?£  fl'i*  ^  " 

f-£ 

Notice  that  the  result  is  independent  of  Ci*  Since  the 
differential  equation  is  satisfied  inside  of  the  cone 
i  >  n  there  are  no  sources  inside  of  this  Mach  cone  and 
since  the  result  is  independent  of  there  are  no  sources 
at  the  Mach  cone. 


The  expression 


■'if'  * 

gives  the  total  outflow  of  mass  due  to  the  perturbation 
through  the  portion  of  the  surface  lying  between  n  =  -e  and 
n  =  ri".  The  total  outflow  of  mass  up  tends  to  infinity  as 
rf  >  0.  (For  the  two  dimensional  case  one  has  a  concentrated 
finite  additional  outflow  along  the  corresponding  characteristic.) 
For  larger  values  of  n  the  total  additional  outflow  becomes 
smaller,  at  £  =  it  reaches  the  value  1.  The  perturbation 
of  the  mass  flow  vector  inside  the  cone  therefore  causes 
inside  a  reduction  of  the  outflow.  The  expression  (A40)  is 
the  fundamental  solution  for  the  three  dimensional  problem. 

Its  form  is,  of  course,  too  singular  to  serve  as  the 
basis  for  a  numerical  procedure.  The  singularity  is  brought 
into  a  less  severe  form  if  one  considers  source  distributions 
and  carries  out  some  of  the  necessary  integrations  beforehand. 


Consider  a  three  dimensional  problem  and  assume  that 
the  source  strength  is  constant  along  the  z  axis.  In  this 
manner  one  will  obtain  the  two  dimensional  potential.  Using 
Eg.  (A40)  one  computes 


■»/*  /5$v**  "/■?  / 


//V-e"; 


'4 


/ft 


t*-w 

This  is  indeed  the  result  Eq.  (A35) . 

To  see  what  happens  if  the  source  density  in  the  y 
direction  is  not  constant,  we  consider  an  example  with  a 
periodic  source  density  in  the  z  direction.  Consider  first 
the  elliptic  problem 

fx*  +  0yy  +  4 it  m  0 

Periodicity  of  the  source  density  in  the  y  direction  implies 
periodicity  of  the  potential.  Accordingly,  we  set 


One  obtains 


*■  Mfoft  > 

2  +  2  -  -  o 

T*x  "yj  T 

and  in  polar  coordinates 

fry  v  rfr  +  ?*£ 


6& 


Now  we  set 

fifr-j  &)  *  fa) 

Then  one  obtains 

£ty  +v2r- 

This  equation  is  solved  by  a  linear  combination  of  H^livr) 

(2)  ° 
and  H  (ivr)  .  For  r  -+  °°  the  solution  must  tend  to  zero; 
o 

obviously  the  effect  of  sources  periodically  distributed 
along  the  z  axis  will  cancel  at  a  large  distance  from  the 
z  axis.  Hence 

J  •  Const  i>r) 

The  constant  is  determined  by  the  local  source  density  at 
r  =  0. 


has 


The  supersonic  problem  is  treated  analogously.  One 

7  >  ..r 

-  Ja*  ~»r  -  0 


2  2  2  ~  _ 
and  with  r  =  x  -  y  and  <J>  =  ij>(r) 


fry  +  0 
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This  equation  is  solved  by 

4  =  J0  r*/^) 

Consider  this  expression  at  a  fixed  x.  The  argument  of  JQ 
varies  between  0  for  the  Mach  cone  and  x  for  the  x  axis. 

The  variation  is  small  if  x  is  small.  The  jump  of  <f> 
at  the  characteristic  r  =  0  is  the  same  for  all  values  of 
x;  namely  JQ(0).  For  large  values  of  x,  the  expression  is 
an  oscillatory  function  of  y.  This  result  has  obvious 
implications  for  the  choice  of  a  grid  in  three-dimensional 
problems . 


Figure  7.  Weight  Functions  for  Case  3b.  Figure  7a  is  also 
the  Weight  Function  for  Cases  5b  and  6b. 


Figure  10.  Approximation  for  $  obtained  for  v  =  0  by  <J>.  =  0, 
4)^  *  *  1  and  exact  function  <J>  for  y  ■  2ir. 
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Ratio  of  Approximate  and  Exact  Wave  Velocities 
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f  Approximate  and  Exact  Wave  Velocities 


io  of  Approximate  and  Exact  Wave  Velocities  Comparison 
Formulae  Based  on  Third  Degree  Splines. 


Ratio  of  Approximate  and  Exact  Wave  Velocities  ,  Comparison 
Between  Second  Degree  and  Third  Degree  Shape  Functions. 


Figure  24.  Amplification  Factor  p  in  the  Complex  p  Plane 

and  Corresponding  Values  of  exp(iyh.),  for  Linear 
Shape  Functions  and  Constant  Weight ^Between 
tk-e  and  tm-e. 
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gure  26.  Amplification  Factor  p  in  the  Complex  Plane  for 
Different  Values  of  =  vh^.  The  arrows  give 

the  direction  of  increasing  vht.  In  part  of  the 
figure  the  points  p  are  connected  with  the  ideal 
values  exp(ivht) .  Third  degree  shape  functions, 
weight  functions  constant  for  0  <  ch*  and 

chfc  <  At  <  hfc,  c  =  3/4.  * 
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Figure  27 .  Amplification  Factor  p  in  the  Complex  Plane  for 
Different  Values  of  A 1/2  =  vht.  The  arrows  give 
the  direction  of  increasing  vhf  In  part  of  the 
figure  the  points  p  are  connected  with  the  ideal 
values  exp(ivht) •  Third  degree  shape  functions, 
weight  functions  constant  for  0  <  At  <  ch^  and 
cht  <  At  <  h^.,  c  *  1. 
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of  Approximate  and  Exact  Wave  Velocities  Third  Degree  Shape 
ons.  Weight  Functions  Give  Exact  Results  for  Very  Long'  Waves 
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Figure  30.  Characteristic  Coordinates  and  Corrections  to  the 
Mass  Flow  Vector  for  Two-Dimensional  Fundamental 
Solutions  in  Linearized  flow  with  M  =  /?. 


Figure  31.  Characteristic  Coordinates  and  region  R  in  the 
CiH-plane  for  an  Axi symmetric  Linearized  Flow 
\fith  M  *  /T. 


